Division by a wrong value in an algorithm?

[max-privato]

Description

I've opened the following ticket: https://github.com/PowerGrids/PowerGrids/issues/52 where the following error occurs: IEEE14busShort4_12jac.c:14213: Invalid root: (0)^(-0.5).

After that I realized that the possible cause of the error is that OM uses a wrong value for R in the denominator of the third row of this algorithm of PowerGrids BusFault:

algorithm
   when pre(fault) then
     Y := 1/Complex(R, X);
   end when;
   when not pre(fault) then
     Y := Complex(0);
   end when;

if this is true, the problem is not a PowerGrids issue, but a OM's. Note that if I add Modelica.Utilities.Streams.print("*******\n" + "R=" + String(R) + ", X=" + String(X)); just before Y := 1/Complex(R, X); the output shows a correct, non-zero value for R. So, the hypothesis is that somehow the C code has an internal flaw that causes an incorrect division by zero.

Unfortunately, I cannot check with Dymola because this model, as well as the original PowerGrids.Examples.IEEE14bus.IEEE14busStaticNetwork, does not initialize in Dymola.

Steps to Reproduce

see on https://github.com/PowerGrids/PowerGrids/issues/52

Expected Behavior

this division by sqrt(0) should not occur

Screenshots

Version and OS

• OpenModelica Version: [e.g. omc --version or Help->About OMEdit from OMEdit]
• OS: [e.g. Windows 10, 64 bit]
• Versions of used Modelica libraries if applicable

Additional Context

[phannebohm]

The division by zero occurs while evaluating some Jacobian entry for Load4.port.U, resulting in

Load4.port.U.$pDERNLSJac66.dummyVarNLSJac66 = 0.5 * (Load4.port.u.re ^ 2.0 + Load4.port.u.im ^ 2.0) ^ (-0.5) * 2.0 * (Load4.port.u.re * Load4.port.u.re.$pDERNLSJac66.dummyVarNLSJac66 + Load4.port.u.im * Load4.port.u.im.$pDERNLSJac66.dummyVarNLSJac66)

At first glance this looks correct to me since the model says

  Real Load4.port.U = Modelica.ComplexMath.'abs'(Load4.port.u);

and ${d\over dt}\Vert u\Vert = (u \cdot \dot{u} )\Vert u\Vert^{-1}$ (using real vector space notation), so we're differentiating correctly, yay!

This suggests that Load4.port.u = 0 + 0i and indeed that's what happens, starting at time=1.2 which is exactly when fault switches from true to false. This leads to 0 = busShort.y = busShort.i = busShort.port.i = busShort.port.S = busShort.port.P = ... I didn't track it all down but lots of (Complex and Real) variables become zero, just from binding equations and pure products. Maybe that gives you a hint? Anyways, that algorithm block seems to not be the source of the problem.

[max-privato]

Thank you for the analysis.

I hoped that the algorithm were the source of the problem since it would ease finding the solution! From the electric point of view, when fault becomes false, the busFault current "i" becomes zero, from

Y:=Complex(0);
and
i=Y*v;

This gives no physical reason for the voltage to become zero. Indeed, busFault is used in other PowerGrids.ENTSOE.TestCase3 voltage recovers immediately after fault removal, as it should also here.

So there is some reason, currently unknown, that causes the bus voltage to become zero when fault becomes false, which then causes the division by zero we see. This is really strange, since fault removal just changes the value of Y, and does nothing else.

[casella]

I think @phannebohm's analysis is correct. In fact, the problem arises after the fault has cleared, when unfortunately the nonlinear solver latches on a bogus solution with extremely low voltage and high current, instead of getting back to normal operation. During nonlinear iterations, it can probably happen that the voltage actually gets to zero, causing the problem with Jacobians involving the modulus of complex power, which involve raising the voltage to the power of -0.5.

I believe this is a modelling problem, I think we can close this ticket and continue the discussion in PowerGrids/PowerGrids#52