Powerflow convergence

[WinfriedL]

I have several test grids where the powerflow does not converge out of the box, while Integral converges without any sign of difficulties. While I cannot rule out the possibility that I have an error in porting the grid, usually I get near perfect agreement between Integral and pandapower (if it converges).

I implemented a few different additions to the pandapower code to help me in the non-convergent cases which I would like to discuss.

1. I implemented a very simple damped newton-raphson by multiplying the correction term with a constant factor smaller than 1, see e.g. https://github.com/WinfriedL/pandapower/blob/23819f7080245a965fe7d9077d81775532b6dced/pandapower/pf/newtonpf.py#L114-L115 This helps convergence in some of the more problematic cases, although it slows the convergence rate. In principle this could be improved by only reducing the factor when needed, e.g. when voltage magnitudes are unreasonable high. At least one could implement an option to set this damping factor, at the moment I have it hardcoded in newtonpf.py. Would this (or an improved version) be in any way interesting for upstream pandapower?

2. Usually I use a flat start for the powerflow calculations, because in our cases it seems to be better than a dc start. However, since barely any bus voltage magnitude is at 1 pu in the results, an increased voltage magnitude for the flat start to about 1.07 to 1.1 pu helps the convergence in our cases. At the moment there are two possible ways to do this in pandapower. I "faked" an res_bus dataframe with my own flat start values and did a powerflow with init='results', and I hardcoded a higher VM in build_bus.py. Both ways are possible, but not very straightforward. Maybe one could add two (optional) columns to the bus dataframe with initial values for the voltage magnitude and angle?

3. To test the above additions and help me analyse convergence issues, I keep all voltage results from the iterations and store them in ppci["internal"]["Vm_it"] and ppci["internal"]["Va_it"] to plot them later on, see e.g. https://github.com/WinfriedL/pandapower/blob/23819f7080245a965fe7d9077d81775532b6dced/pandapower/pf/newtonpf.py#L126-L127 Of course this increases memory consumption and slows down execution (have not tested if it is in any way significant), but it might be useful for others as well?

4. In one particular case I had the problem that depending on the used parameters for damping and initial voltages I got two different powerflow results. Both were (according to pandapower) fully converged, but the voltages (magnitude and angle) differed significantly.
   First I came across this with the damping method, where I got the right result (compared to Integral) with a damping of about 0.5, and a different result with a damping of 0.3. But I also see the same results with different flat start values and no damping, where I have ranges of values with one of the results and ranges of values with the other. Do you have any idea where these two different seemingly correct results come from? In principle there should be two mathematical correct solutions, is this a possible explanation? Is there any way I can check this? Did you ever come across such a behaviour?

[lthurner]

1. There is an algorithm called Levenberg-Marquardt that works with a damping factor that is adjusted in every step and is supposed to work well for ill-conditioned systems. I am not familiar with the details of how the algorithm works, but that might be something worth looking into. GridCal has an implementation of the algorithm: https://github.com/SanPen/GridCal/blob/master/UnderDevelopment/GridCal/Engine/Numerical/JacobianBased.py#L280. Since it is also based on the pypower data structure, it might be relatively straightforward to adopt this code.

2. I am in favor of giving an option in runpp that allows to use a different value than 1 for flat start (e.g vm_pu_flat=1.). If one wants to define a different starting value for each bus, I think "faking" a res_bus dataframe is an OK solution, and its not worth having a specific initialization mode and an extra column in bus column to achieve the same functionality. Maybe a convenience function that adds the bus voltages to the res_bus table could be added the toolbox?

3. This seems useful for debugging but nothing that should be done by default. We could discuss if it is worth to include an option in runpp for this...

4. I know the power flow problem can theoretically have multiple solutions, but working with mainly radial systems I have never personally had this problem. Maybe someone else knows more about that.

[FlorianShepherd]

From your comments I understand that speed is not the main issue, but convergence is crucial. In my opinion an implementation of a HELM solver would be an elegant solution. It guarantees to find the correct solution, but takes some more time to convergence in comparison to NR.

Here are links to two papers about the HELM Solver: https://ieeexplore.ieee.org/document/6344759/ https://ieeexplore.ieee.org/document/7352383/

Also GridCal has a implementation in progress: https://github.com/SanPen/GridCal/wiki/The-HELM-Algorithm

I tried copying it, but it didn't work (at least in half an hour). I'm not sure how to get Yseries and if Yshunt is the same as in pypower.

[WinfriedL]

Thanks for your replies!

There is an algorithm called Levenberg-Marquardt that works with a damping factor that is adjusted in every step and is supposed to work well for ill-conditioned systems.

There are indeed several algorithms which make use of a damping factor. I just implemented something with minimal code changes by just using a constant factor.

I am in favor of giving an option in runpp that allows to use a different value than 1 for flat start (e.g vm_pu_flat=1.)

I think this might be helpful. On the other hand, if the toolbox had some helper functions to create the res_bus table, it might not be needed. The problem with 'faking' the res_bus table was for me to set all voltage controlled buses to their respective voltage setpoint. I don't know (and have not checked) if it is a problem for pandapower if the initial voltage value differs from the generator voltage? In addition, I tried to set all bus voltages with extended wards to the extended ward voltage, because in our cases the impedance between the internal end external bus is usually low and I assumed it would be a better starting point. For both cases, generators and extended wards, and maybe for other use cases as well, you would have to know which pandapower buses are fused together by switches, because all those buses should have the same initial voltage. I used the bus lookup table for that, but in general you don't have acces to the lookup until you run the powerflow.

This seems useful for debugging but nothing that should be done by default.

I agree, I used it a lot for debugging. Just to illustrate, here are two voltage results with different flat start values, the top one uses vm0=1.06, the bottom one vm0=1.07. As you can see, both are looking converged, but the results differ.

I know the power flow problem can theoretically have multiple solutions, but working with mainly radial systems I have never personally had this problem.

Me neither. It worries me that I have two different 'converged' solutions for the same testcase, depending on initial voltages (which seem all reasonable) or damping.

From your comments I understand that speed is not the main issue, but convergence is crucial.

Well yes. If it does not converge, it does not help me if it is fast ;) My two main concerns are getting convergence at all (while Integral converges without problems, and does not see any instability issues), and don't getting two different results

[lthurner]

I think this might be helpful. On the other hand, if the toolbox had some helper functions to create the res_bus table, it might not be needed.

Since the DC initialization only calculates voltage angles, it would make sense to be able to combine a starting value for the voltage magnitude with a DC initialization. That would not be possible with init_results (except by carrying out a DC power flow first and writing the angle results in the result table, which is not exactly straightforward). So I think a new parameter that allows to define a starting voltage magnitude that is used for flat and DC initialization would still make sense.

The problem with 'faking' the res_bus table was for me to set all voltage controlled buses to their respective voltage setpoint. I don't know (and have not checked) if it is a problem for pandapower if the initial voltage value differs from the generator voltage?

Yes all PV/slack buses have to be initialized with their voltage setpoint in order for the power flow to converge, because the newton-raphson never changes those voltages, so it would be impossible to reach the setpoint at those buses. For flat start, that is handled internally, but for init_results I am not sure? The voltage magnitudes at PV/slack buses should be overwritten with the setpoint valuse even when init_results is used, because if there is a timeseries simulation that initializes with the result of the previous power flow, but the voltage setpoint of a generator/slack/ward changes between timesteps, that might lead to incorrect results or non-convergence. I am not sure how this is currently implemented, but we should definitely check that...

In addition, I tried to set all bus voltages with extended wards to the extended ward voltage, because in our cases the impedance between the internal end external bus is usually low and I assumed it would be a better starting point.

That seems like generally a good assumption, I think its worth considering doing this by default for extended wards.

For both cases, generators and extended wards, and maybe for other use cases as well, you would have to know which pandapower buses are fused together by switches, because all those buses should have the same initial voltage. I used the bus lookup table for that, but in general you don't have acces to the lookup until you run the powerflow.

Instead of using the lookup, you can also solve this problem with a graph search, which also allows you to find which buses are galvanically connected. In fact I just remembered that I already implemented a function that uses graph searches to figure out starting point voltage magnitudes and angles by graph searches: https://github.com/lthurner/pandapower/blob/develop/pandapower/topology/graph_searches.py#L454. That function assumes a radial grid, and just starts from each external grid and initializes all connected buses to the external grid set point. It also considers the ratio and angle shift of transformers along the way. For meshed grids, one would have to think of a different algorithm, but the basic idea could also be applied.

[lthurner]

Regarding the voltage setpoint at slack/PV buses, I just checked and the setpoint values are always written in the ppc["bus"] matrix in the build_gen functions, for example here for the slack buses: https://github.com/lthurner/pandapower/blob/develop/pandapower/build_gen.py#L235. So regardless of how the voltages are initialized, the PV/slack buses are always initialized with the setpoint (which is the correct behaviour imho), and you should not have to worry about those values when "faking" a res_bus for initialization.

[WinfriedL]

You are right, I just tried. Don't know why, but during my early testing I had problems when setting all voltage magnitudes in res_bus to a constant value. Now it works as expected. This makes creation of a (constant) flat start a lot easier. Only in case one wants to sepcify individual start voltages one would have to deal with switches. I have no idea what pandapower does when several fused buses have different initial voltages... without looking at the code I would expect that only one of the values gets used...

by carrying out a DC power flow first and writing the angle results in the result table, which is not exactly straightforward

Maybe not straightforward, but relativly easy to test...

[WinfriedL]

I have several updates...

1. While I found the idea to initialise the angles with a dc powerflow and use an 'improved' flat start for the magnitudes interesting, it does not really change anything. The convergence behaviour (all the mentioned symptoms) is basically the same between a dc start and va=0. The convergence rate is also the same.

2. I could confirm that in my testcase I run into the multiple possible solutions for the powerflow equations, so from a mathematical/numerical standpoint there is nothing wrong.

3. I could reproduce the effects I saw with pandapower in Integral by overriding Integrals (unknown) default flat start. If I use a flat start with vm=1 and va=0, both pandapower and Integral do not converge. If I use a flat start with vm=1.06, both pandapower and Integral do converge, but get the technically wrong result. And with a flat start with vm=1.07, both converge and give the same (right) results. So there is no problem with either my port or pandapower.

4. This leaves me with the problem to find a) the best flat start for our use cases and b) a way to automatically identify unreasonable solutions and throw a warning that case. Any ideas are welcome.

5. You/we can still discuss which of my code additions might be interesting to implement.

[jhmenke]

Usually I use a flat start for the powerflow calculations, because in our cases it seems to be better than a dc start. However, since barely any bus voltage magnitude is at 1 pu in the results, an increased voltage magnitude for the flat start to about 1.07 to 1.1 pu helps the convergence in our cases. At the moment there are two possible ways to do this in pandapower. I "faked" an res_bus dataframe with my own flat start values and did a powerflow with init='results', and I hardcoded a higher VM in build_bus.py. Both ways are possible, but not very straightforward. Maybe one could add two (optional) columns to the bus dataframe with initial values for the voltage magnitude and angle?

Maybe the initialization option init="slack" from the state estimation can be ported to the power flow calculation (should be simple). It detects initial voltage angle changes due to transformers and sets initial bus voltage magnitudes to the value set at the connected slack bus.

Maybe you can try it quickly by setting res_bus to the result of https://github.com/lthurner/pandapower/blob/f14386ceccd21ae24a2e51962767aced7a8fd80e/pandapower/topology/graph_searches.py#L454 ?

[WinfriedL]

Maybe the initialization option init="slack" from the state estimation can be ported to the power flow calculation (should be simple). It detects initial voltage angle changes due to transformers and sets initial bus voltage magnitudes to the value set at the connected slack bus.

Thanks for the suggestion, but it won't really help in our networks.

• In our large meshed networks, it would basically set all voltages to the setpoint of the slack bus. Best case would be a better flat start than vm_pu=1.0, but still constant.
• We have multiple voltage controlled generators with different voltage setpoints in the same meshed network. Using only the slack voltage for all other buses might lead to even worse convergence.
• The slack bus in our meshed network has no real technical meaning appart from being the reference bus with va_degree=0. In fact we choose the slack bus from all PV-buses based on the parameters of the connected generators, so it would not be well defined in my opinion. The flat start would depend on which bus gets chosen as the slack and its voltage setpoint. In this particular testcase the slack bus has a voltage setpoint of vm_pu=1.0.
• I guess it does make much more sense in radial networks where the slack bus (external grid) has a real meaning.

[jhmenke]

Okay, that makes sense. We almost exclusively deal with distribution grids.

This leaves me with the problem to find a) the best flat start for our use cases and b) a way to automatically identify unreasonable solutions and throw a warning that case. Any ideas are welcome.

How do you evaluate that the solution found by using flat 1.06 p.u. is not valid but the 1.07 p.u. solution is valid? If there are objective, measurable criteria then a simple for loop with different flat start positions and an automatic evaluation afterwards could be a starting point.

[WinfriedL]

How do you evaluate that the solution found by using flat 1.06 p.u. is not valid but the 1.07 p.u. solution is valid?

• I compared to Integral and the solution with 1.07 p.u. in in perfect agreement, while the solution with 1.06 is not
• Several buses have lower resulting voltages in the invalid solution
• Transformers connected to these buses have huge reactive power losses, (10 GVar)
• That already hints at the invalid solution with lower voltages and higher currents
• Finally, I calculated PV-curves for some of these buses and the results were obvious:

Problem is: I don't see any of this as a way to determine automatically (and in reasonable time) if the solution is valid or not. Appart from the reactive power losses, but I don't know if that happens necessarily in the invalid solution.

[WinfriedL]

I think I found a way to get reasonable way to get initial flat start values for my cases, based on the generators in the network. Huge thanks for all the topological functionality in pandapower!

So at the moment I don't need a different or modified algrotihm. I would suggest closing this issue and maybe open new individual issues for any of the improvements you are interested in, like

• keeping results per iteration for debugging
• finding/adding some kind of stability test for the solution
• adding additional algorithms like Levenberg-Marquardt or HELM

Like I said, I don't have any urgent need for any of this, but I definitly would help testing. As for the first one, I could make a pull request if you are interested.

[lthurner]

Great that you found a way to make it work with graph searches. How did you solve the problem now? I imagine this might also be relevant for other users, maybe you could provide that function as a helper function for the toolbox?

If you could make pull request for the voltage result, that would be appreciated. The logging should be disabled by default so that it doesn't affect performance, so we would probably need a new parameter in runpp?

I am going to open new issues for other algorithms

[WinfriedL]

If you could make pull request for the voltage result, that would be appreciated. The logging should be disabled by default so that it doesn't affect performance, so we would probably need a new parameter in runpp?

I will make a pull request, but probably not until sometime next week. I will test if it affects performance as well. If there is no measurable influence, maybe it is not necessary to add a new parameter.

[WinfriedL]

I have done some tests on the performance impact:

NR steps	runtime increase [%]	needed memory [MB]
5	-1.1	0.46
9	1.2	0.78
11	1.4	0.92
45	8.9	3.6

The result for 5 NR steps is probably just statistics, the needed memory is the size of the two numpy matrices together.

So for single calculations with normal number of NR steps the impact is small. But even a 1-2 % increase in calculation time could justify an additional parameter, if one has to do lots of calculations...