Questionable tearing in PowerSysPro leads to initialization failure

[Francesco Casella]

Please check ​PowerSysPro.Examples.OneProdLoadOneLineOneLineOneLoad. The initialization fails because of a division by zero in a torn equation.
The zero term is line2.Y.re, which depends on line2.G = 0.

This choice of parameter is perfectly legitimate, the line is purely resistive and has zero shunt admittance, which should lead to a solvable system of equations. In fact, by setting --tearingMethod=noTearing the system can be initialized successfully.

I guess we should improve the tearing algorithm to avoid this kind of situation. This may actually have a positive influence on other failing models in the testsuite.

[Karim Adbdelhak]

You can prevent division by parameters with --tearingStrictness=veryStrict. We can make that option default if you like. Seems like it should solve this problem.

[Francesco Casella]

Yeah, that's one option.

Maybe we should explore if we can figure out when the division by zero only happens if the system is singular (in which case it is fine to keep it) or not (in which case it would be better to avoid it). If we just check for structural singularity, we should keep it if the element in question of the Jacobian is the only non-zero one in that column and row.

[Francesco Casella]

This issue is currently being worked on in ​PR 6969. As soon as that PR is merged, I expect most, if not all, of PowerSysPro models to work.

@jean-philippe, at that point, you will need the 1.17.0-dev nightly build to run your library. We may release a 1.17.0-dev.01 stable version so you have a stable release to run it until the release of 1.17.0 in Feb 2021.

[jean-philippe.tavella@…]

@Francesco: thank you for the information. According to our internal plans at EDF, I think better to wait till the delivery of version 1.17.0 in Feb. 2021.

In the meanwhile with 1.16.0, do we preferentially have to use:
--tearingStrictness=veryStrict
or --tearingMethod=noTearing

[jean-philippe.tavella@…]

Replying to casella:

This issue is currently being worked on in ​PR 6969. As soon as that PR is merged, I expect most, if not all, of PowerSysPro models to work.

@jean-philippe, at that point, you will need the 1.17.0-dev nightly build to run your library. We may release a 1.17.0-dev.01 stable version so you have a stable release to run it until the release of 1.17.0 in Feb 2021.

I tested 1.17.0-dev-315 and I saw the flags OMCFlags="-d=initialization --TearingMethod=minimalTearing".
But when I reinit the options, these flags disappeared. Is it intended?

[Francesco Casella]

Replying to jean-philippe.tavella@…:

@Francesco: thank you for the information. According to our internal plans at EDF, I think better to wait till the delivery of version 1.17.0 in Feb. 2021.

OK. Sorry for the delay in the reply.

In the meanwhile with 1.16.0, do we preferentially have to use:
--tearingStrictness=veryStrict
or --tearingMethod=noTearing

In principle veryStrict should be enough to make sure that no division by zero is ever obtained, though I'm not 100% sure. Maybe @Karim or @AnHeuermann can comment on that.

However, I am convinced that as long as we can now use sparse solvers, there is really no big advantage at using tearing for linear system in terms of performance. So, our plan is to avoid tearing for linear systems, and rely on ​-lssMaxDensity and ​-lssMinSize to automatically switch to the sparse solvers. Maybe we should revise those limits, the problem is, we don't have a good test suite to base the decision on actual performance data.

That said, until ​PR 6969 is merged in , you should use minimalTearing, that still applies tearing for discrete unknowns. If you set noTearing, systems involving discrete unknowns will no longer be solvable.

[Francesco Casella]

Replying to jean-philippe.tavella@…:

I tested 1.17.0-dev-315 and I saw the flags OMCFlags="-d=initialization --TearingMethod=minimalTearing".

I guess you see them because they are set in custom __OpenModelica_commandLineOptions and __OpenModelica_simulationFlags annotations in the model. That's the safest way to set flags, because it is embedded in the Modelica source code and it works on every machine and installation out of the box.

You can see some examples of that in the PowerGrids.Examples package.

But when I reinit the options, these flags disappeared. Is it intended?

What do you mean by "reinit the options" exactly?

[jean-philippe.tavella@…]

Replying to casella:

Replying to jean-philippe.tavella@…:

I tested 1.17.0-dev-315 and I saw the flags OMCFlags="-d=initialization --TearingMethod=minimalTearing".

I guess you see them because they are set in custom __OpenModelica_commandLineOptions and __OpenModelica_simulationFlags annotations in the model. That's the safest way to set flags, because it is embedded in the Modelica source code and it works on every machine and installation out of the box.

You can see some examples of that in the PowerGrids.Examples package.

But when I reinit the options, these flags disappeared. Is it intended?

What do you mean by "reinit the options" exactly?

I work under Windows, via the GUI OMEdit. Reinit is possible with the command Tools->Options->Simulation and then the button 'Réinitialiser' in French (probably Reinit).
Then it is proposed to save an omedit.ini file with all the flags saved.
An idea on the way to reactivate this file, or another one if wanted by user?

[jean-philippe.tavella@…]

That said, until ​PR 6969 is merged in , you should use minimalTearing, that still applies tearing for discrete unknowns. If you set noTearing, systems involving discrete unknowns will no longer be solvable.

I will use --TearingMethod=minimalTearing till PR6969 is merged.

[Francesco Casella]

Replying to jean-philippe.tavella@…:

I work under Windows, via the GUI OMEdit. Reinit is possible with the command Tools->Options->Simulation and then the button 'Réinitialiser' in French (probably Reinit).

I guess it is "Reset" in English. This means that all settings are reset to factory standard.

Then it is proposed to save an omedit.ini file with all the flags saved.

If you do that, you are overwriting the settings that you may have put in the model earlier, e.g. to set tearingMethod

An idea on the way to reactivate this file, or another one if wanted by user?

Once you've overwritten and saved, they are gone. Of course if you use a revision control system you can always roll back to the previous version of the file.

[Francesco Casella]

Rescheduled to 1.18.0

[Francesco Casella]

Following our discussion at the devmeeting today, I double-checked the ​documentation of strict tearing, and it says

(Robust tearing rules by consideration of the partial derivative. Allows to divide by parameters that are not equal to or close to zero.)

According to the documentation, division by such parameters should be avoided when strict tearing (the default) is performed. Apparently something goes wrong. @Karim, can you please double-check this?

[Karim Adbdelhak]

Everything seems to work for me using OpenModelica 1.17.0 and OMCompiler v1.18.0-dev.204+g1a8d292348. Following simulation result and dump of the torn component for initialization seem fine for me. There is an initial assert for some value, but that seems to be only triggered during the search for a viable initialization.

Maybe this has already been fixed? Can someone double check my results?

record SimulationResult
    resultFile = "/home/kab/Repo/FH/promotion/modelicaLibs/powersyspro/test/PowerSysPro.Examples.OneProdLoadOneLineOneLineOneLoad_res.mat",
    simulationOptions = "startTime = 0.0, stopTime = 1.0, numberOfIntervals = 500, tolerance = 1e-06, method = 'dassl', fileNamePrefix = 'PowerSysPro.Examples.OneProdLoadOneLineOneLineOneLoad', options = '', outputFormat = 'mat', variableFilter = '.*', cflags = '', simflags = ''",
    messages = "assert            | warning | The following assertion has been violated during initialization at time 0.000000
|                 | |       | (-pv.P) > (-pv.Pmax)
assert            | warning | >>> Active power is below the maximum acceptable value for PowerSysPro.Examples.OneProdLoadOneLineOneLineOneLoad.pv
LOG_SUCCESS       | info    | The initialization finished successfully without homotopy method.
LOG_SUCCESS       | info    | The simulation finished successfully.
",
    timeFrontend = 0.165857223,
    timeBackend = 0.079235921,
    timeSimCode = 0.007037741,
    timeTemplates = 0.010965698,
    timeCompile = 1.149814192,
    timeSimulation = 0.040691202,
    timeTotal = 1.454352817
end SimulationResult;
""

################################################################################
 dumpLoops: SORTED COMPONENT 
################################################################################

torn nonlinear Equationsystem:

internal vars (4)
1: pv.P:VARIABLE(unit = "W" )  "Active power flowing the node" type: Real
2: pv.Q:VARIABLE(unit = "var" )  "Reactive power flowing the node" type: Real
3: line2.terminalA.i.re:VARIABLE(flow=true start = 0.0 unit = "A" )  "Real part of complex number" type: Real
4: line2.terminalA.i.im:VARIABLE(flow=true start = 0.0 unit = "A" )  "Imaginary part of complex number" type: Real

residual vars (10)
1: line1.terminalA.v.re:VARIABLE(flow=false start = line2.UNom / 1.732050807568877 unit = "V" )  "Real part of complex number" type: Real
2: line1.terminalA.v.im:VARIABLE(flow=false start = 0.0 unit = "V" )  "Imaginary part of complex number" type: Real
3: line2.terminalB.i.im:VARIABLE(flow=true start = 0.0 unit = "A" )  "Imaginary part of complex number" type: Real
4: line2.terminalB.i.re:VARIABLE(flow=true start = 0.0 unit = "A" )  "Real part of complex number" type: Real
5: line1.terminalB.i.im:VARIABLE(flow=true start = 0.0 unit = "A" )  "Imaginary part of complex number" type: Real
6: line1.terminalB.i.re:VARIABLE(flow=true start = 0.0 unit = "A" )  "Real part of complex number" type: Real
7: line2.terminalB.v.im:VARIABLE(flow=false start = 0.0 unit = "V" )  "Imaginary part of complex number" type: Real
8: line2.terminalB.v.re:VARIABLE(flow=false start = load.UNom / 1.732050807568877 unit = "V" )  "Real part of complex number" type: Real
9: pv.terminal.v.re:VARIABLE(flow=false start = line1.UNom / 1.732050807568877 unit = "V" )  "Real part of complex number" type: Real
10: pv.terminal.v.im:VARIABLE(flow=false start = 0.0 unit = "V" )  "Imaginary part of complex number" type: Real

internal equations (4)
1/1 (1): 3.0 * ((-pv.terminal.v.re) * line1.terminalB.i.re - pv.terminal.v.im * line1.terminalB.i.im) = pv.P   [dynamic |0|0|0|0|]
2/2 (1): 3.0 * (-pv.terminal.v.im) * line1.terminalB.i.re + (-3.0) * (-pv.terminal.v.re) * line1.terminalB.i.im = pv.Q   [dynamic |0|0|0|0|]
3/3 (1): line2.terminalA.i.re + line2.terminalB.i.re = line2.Y.re * (line2.terminalB.v.re + line1.terminalA.v.re) - line2.Y.im * (line2.terminalB.v.im + line1.terminalA.v.im)   [dynamic |0|0|0|0|]
4/4 (1): line2.terminalA.i.im + line2.terminalB.i.im = line2.Y.re * (line2.terminalB.v.im + line1.terminalA.v.im) + line2.Y.im * (line2.terminalB.v.re + line1.terminalA.v.re)   [dynamic |0|0|0|0|]

residual equations (10)
1/1 (1): pv.terminal.v.re = 0.5773502691896258 * pv.UNom / sqrt(1.0 + (pv.Q / pv.P) ^ 2.0)   [dynamic |0|0|0|0|]
2/2 (1): pv.terminal.v.im = 0.5773502691896258 * pv.UNom * pv.Q / (sqrt(1.0 + (pv.Q / pv.P) ^ 2.0) * pv.P)   [dynamic |0|0|0|0|]
3/3 (1): line1.terminalB.i.re - line2.terminalA.i.re = line1.Y.re * (pv.terminal.v.re + line1.terminalA.v.re) - line1.Y.im * (pv.terminal.v.im + line1.terminalA.v.im)   [dynamic |0|0|0|0|]
4/4 (1): line1.terminalB.i.im - line2.terminalA.i.im = line1.Y.re * (pv.terminal.v.im + line1.terminalA.v.im) + line1.Y.im * (pv.terminal.v.re + line1.terminalA.v.re)   [dynamic |0|0|0|0|]
5/5 (1): line1.terminalA.v.im - pv.terminal.v.im = line1.Z.re * ((-line1.Y.im) * line1.terminalA.v.re - line2.terminalA.i.im - line1.Y.re * line1.terminalA.v.im) + line1.Z.im * (line1.Y.im * line1.terminalA.v.im + (-line1.Y.re) * line1.terminalA.v.re - line2.terminalA.i.re)   [dynamic |0|0|0|0|]
6/6 (1): 3.0 * (-line2.terminalB.v.im) * line2.terminalB.i.re + (-3.0) * (-line2.terminalB.v.re) * line2.terminalB.i.im = load.Q   [dynamic |0|0|0|0|]
7/7 (1): 3.0 * ((-line2.terminalB.v.re) * line2.terminalB.i.re - line2.terminalB.v.im * line2.terminalB.i.im) = load.P   [dynamic |0|0|0|0|]
8/8 (1): line1.terminalA.v.im - line2.terminalB.v.im = line2.Z.re * (line2.terminalA.i.im + (-line2.Y.re) * line1.terminalA.v.im - line2.Y.im * line1.terminalA.v.re) + line2.Z.im * (line2.terminalA.i.re + line2.Y.im * line1.terminalA.v.im - line2.Y.re * line1.terminalA.v.re)   [dynamic |0|0|0|0|]
9/9 (1): line1.terminalA.v.re - line2.terminalB.v.re = line2.Z.re * (line2.terminalA.i.re + line2.Y.im * line1.terminalA.v.im - line2.Y.re * line1.terminalA.v.re) - line2.Z.im * (line2.terminalA.i.im + (-line2.Y.re) * line1.terminalA.v.im - line2.Y.im * line1.terminalA.v.re)   [dynamic |0|0|0|0|]
10/10 (1): line1.terminalA.v.re - pv.terminal.v.re = line1.Z.re * (line1.Y.im * line1.terminalA.v.im + (-line1.Y.re) * line1.terminalA.v.re - line2.terminalA.i.re) + line1.Z.im * (line1.Y.re * line1.terminalA.v.im - ((-line1.Y.im) * line1.terminalA.v.re - line2.terminalA.i.im))   [dynamic |0|0|0|0|]