incorrect calculation of ZIP-Loads with voltage angles

[rbolgaryn]

In the following test case:

• loads have const_i_percent=100
• transformer has phase shift != 0 (5*30)
• calculate_voltage_angles = True

the load flow results are false:

• negative values in net.res_load
• wrong voltage magnitudes and degrees in net.res_bus

A simple example to recreate the issue is in the file test.p

test.zip

P. S. see test_zip_loads_results() in test_runpp.py

[lthurner]

I tried it with an external grid with a voltage angle instead of the transformer with angle shift, and the same problem occurs. So it does not seem to be related to the transformer model, but rather to how the I part of the ZIP load is calculated from the complex voltage vector.

I changed the test to a more minimal two bus example without transformer as test case for this problem (see commit f576d6f)

[jkrstulo]

My explanation for this issue is the following (I'll refer to the specific example in order to be more clear):
load is set as p_kw=2 with constant current part 100%...in the power flow calculation this is treated as constant current injection at the corresponding bus. Value of the constant current is calculated for the nominal voltage 1 p.u. which results in this case with -0.002+0.j.
The problem is that in the case of considering voltage angles, nominal voltage is not 1.+0.j p.u. but it has angle 100 degrees, that is -0.17364818+0.98480775j
so in this specific case, constant current at the bus and voltage have 100 degrees phase difference (power factor cos(phi)=-0.17) so the load results with the value p_kw=-0.348496 q_kvar=1.968417

So one solution would be to consider initial voltage angle for calculation of the constant current loads. In this case constant current is constant vector, so the power factor changes also slightly with the voltage. This solution is implemented in GridCal, which I initially used for the reference.

Alternative would be to consider constant current magnitude, while keeping fixed power factor. so, active and reactive power is independently calculated as:

Please let me know if you have some other suggestion?

P.S. and sorry for taking me so long to address this issue

[lthurner]

The problem is that in the case of considering voltage angles, nominal voltage is not 1.+0.j p.u. but it has angle 100 degrees, that is -0.17364818+0.98480775j
so in this specific case, constant current at the bus and voltage have 100 degrees phase difference (power factor cos(phi)=-0.17) so the load results with the value p_kw=-0.348496 q_kvar=1.968417

Remember that the voltage angle is not the angle between current and voltage, but between bus voltage angle and reference angle. The cos(phi) at the bus is only determined by the P/Q ratio of the bus injection, the voltage angle should have no impact on it. That is why the result should not be different when the slack angle is changed... I am not sure why the current calculation changes with the voltage angle...

[jkrstulo]

Remember that the voltage angle is not the angle between current and voltage, but between bus voltage angle and reference angle

of course, don't worry I am not mixing these angles...

The cos(phi) at the bus is only determined by the P/Q ratio of the bus injection, the voltage angle should have no impact on it.

I also fully agree...I was not justifying that this is correct, just trying to explain why it happens: constant current is calculated before the power flow here...the problem is that it is considered as constant phasor with phase angle equal to zero..in the mean time voltage angle is set to 100 and you get load with some crazy power factor instead of power factor = 1 (Q is zero)
so in other words voltage angles are considered but current angles are not and this produces these weird results...

from all this I conclude that we should change zip-loads implementation on the way that current magnitude is kept fixed (not current phase) and in respect to that also power factor is fixed, as you can see in my updated comment beforehand (I added alternative solution with equation, please read again last 2 paragraphs)...

[lthurner]

Yes the second solution sounds exactly right, I think that is it. Thanks for taking the time to look into this Jakov!

[jkrstulo]

yes the 2nd solution will be more suitable we agree, although I think we can't say neither one solution is more correct...

...just I wanted to clarify a bit more the first solution and why I did it like that...
(one of the reasons I implemented it this way is because I followed GridCal as a reference...)
...the idea is: constant current load is defined with active and reactive power at nominal voltage...at that voltage constant current is calculated and that constant current (as a phasor) is fixed until the end of the calculation (the same thing is done with constant impedance).
To be more clear let's see the results for the case without considering voltage angles: load is defined as p_kw=2 (q_kvar=0 so power factor is 1)...however the result is p_kw=1.998801 q_kvar=-0.001201...so the power factor is no more 1 because current was fixed to -0.002+0.j and the resulting voltage phasor is vm_pu=0.999401 va_degree=-0.034425 so obviously voltage and current at the bus 2 are not in phase and that is the reason why power factor is no more 1.
to resume: in this solution current phasor is fixed which results with slight variation of the power factor in the solution (this change of power factor is mostly not significant if voltage angles are close to zero, but in case of phase shifting we get more weird results as is reported in this issue).

in the second solution, we would fix current magnitude and leave the power factor fixed (as defined with p_kw and q_kvar. So the implementation would be significantly different in regard to the first solution: instead of calculating current and impedance and fixing them for the calculation, the load value of the voltage-dependent part will be updated in each iteration (that is Sbus will be recalculated in each iteration) with new voltage values, according to:

and that's it...not any current, or impedance calculation is necessary...

[lthurner]

That is a very nice summary and I understand the difference now.

Just from a theoretical point of view I would have thought calculating the complex current as it is now is correct. But with these specific results, the fixed current magnitude looks more correct to me, because I would not expect to get a result of p_kw=-0.35 q_kvar=1.97 when entering p_kw=2 and q_kvar=0. Commercial tools also seem to implement the fixed current magnitude, so it really seems the way to go.

Since it is a completely different approach (calculating Ibus vs. updating Sbus), will it be much work to implement it that way?

#edit: you mentioned impedances also, does this problem also exist for shunts then?

[jkrstulo]

#edit: you mentioned impedances also, does this problem also exist for shunts then?

no it is not such a problem for impedance, because power is calculated from voltage and impedance, and power factor will depend on impedance, irrespective of the voltage angle. (but i'll still double check it through an example)

Since it is a completely different approach (calculating Ibus vs. updating Sbus), will it be much work to implement it that way?

it shouldn't be a big work I guess...it should be much easier to implement than the way it is implemented now...
but however I suggest the following: I will first bug-fix the current implementation to consider voltage angle when calculating constant-current part - so that current phasor is calculated form the flat start which includes voltage angle (flat start initiated by DC power flow in case of transformer phase shift, or just rotated flat start for external grid angle)...and commit it...this is only to have this solution bug-fixed and passed through git...(who knows if we come back to it in the future)
then I would change the solution completely to fixed current magnitude and leave that as ultimate solution...
if you agree?

[lthurner]

yes sounds good!