Adding sensors changes results

[massimo ceraolo]

The Test1.mo needs PowerGrids to run.
It contains models ScTT and ScTT1. The latter is the same as the former with the addition of a sensor and a block to process its output.

The two models give different results for the variable Gen.port.VPu: 0.808 and 0.956, both constant over the simulation interval.

The system during the simulation interval is linear and time-invariant; it should therefore have just one equilibrium point, corresponding to VPu=0.956.

Tested with v1.17.0-dev-195.

[hanshell <hans.hell@…>]

sorry, your models are not identical.
In ScTT you changed the initial-condition for the GEN to:

localInit = PowerGrids.Types.Choices.LocalInitializationOption.PQ,

which results in an underdetermined system (which is unstable)

Check of Test1.ScTT completed successfully.
Class Test1.ScTT has 328 equation(s) and 332 variable(s).
121 of these are trivial equation(s).

model ScTT1 is OK.

[Francesco Casella]

Replying to ceraolo:

The system during the simulation interval is linear and time-invariant; it should therefore have just one equilibrium point, corresponding to VPu=0.956.

The system is linear during simulation (I'm not 100% sure about time-invariant, due to the Park transformation within the generator), but it's not during initialization. You can check that easily by running the simulation with the declarative debugger (blue bug icon) and by checking how the initial equations are solved.

Please be careful with the initialization options. As discussed in the tutorial and in the Modelica Conference paper there are three initialization modes, leading to different initialization problems, depending on your needs. They need to be applied consistently to the entire system, otherwise meaningless initialization problem may occur.

I would suggest you to set all the start values from the power flow correctly and to see if this solves the problem. If it does, then it is just a Newton method initialization issue, and we can trace back why exactly it fails when you don't set all start values properly. If it doesn't, then we may have some more serious problems in the equation handling, and we'll investigate that.

[massimo ceraolo]

Replying to casella:

Replying to ceraolo:

The system during the simulation interval is linear and time-invariant; it should therefore have just one equilibrium point, corresponding to VPu=0.956.

The system is linear during simulation (I'm not 100% sure about time-invariant, due to the Park transformation within the generator), but it's not during initialization. You can check that easily by running the simulation with the declarative debugger (blue bug icon) and by checking how the initial equations are solved.

The examples I supply here seem really strange because of the following:

1. two systems differing only because the second has one sensor more give different initial values, not only in initialisation, but also during the transient
2. the system modelled is physically linear (in a sense also time-invariant, but better not to discuss this here). It would be non-linear if I wanted to fix active and reactive powers, i.e. setting-up a power flow problem. But here I set only excitation voltage and input mechanical power, while active and reactive terminal powers are just initial guesses, as well as UStart and UPhaseStart.
3. if the system is linear, only one stationary solution exists; but the two systems shown here stay both in steady state (all constant values), but with different values from each other

I know that I used bad values for start values of the machine (UStartm UPhaseStart, PStart, QStart), chosen I don't remember how. I left them there on purpose, and chose not to make a load flow to find better initial values, to show that there is a problem.

In fact, the points 1..3 are valid even if I give bad values for the initialization process: if only an initial point exists, the solver should either find it or not converge; In the unfortunate case it finds a bad initial point, the subsequent transient should not show all constant values.

Well, maybe putting forward an issue in this way is not good for developers, which need simpler MWE's Maybe I can close the ticket for now, until when I find a simpler example to show the issue.

@casella, just one more question: all the initialisation techniques in PowerGrids assume that all inter-connected generators share the system frequency, isn't it?

[Francesco Casella]

Replying to ceraolo:

1. two systems differing only because the second has one sensor more give different initial values, not only in initialisation, but also during the transient

Unfortunately, with equation-based modelling the addition of one component can have unpredictable consequences on the structure of the generated code. If the convergence is critical, this may trigger the solver to converge to another solution. Welcome to the dark side of equation-based modelling...

2. the system modelled is physically linear (in a sense also time-invariant, but better not to discuss this here). It would be non-linear if I wanted to fix active and reactive powers, i.e. setting-up a power flow problem. But here I set only excitation voltage and input mechanical power, while active and reactive terminal powers are just initial guesses, as well as UStart and UPhaseStart.

I get that, but if you initialize in steady-state, the machine angle theta is also obviously an unknown of the initialization problem, and it causes the system to become nonlinear because of Park's transformation. Once you've locked onto one solution at initialization, the subsequent simulation sticks to that by continuity.

In fact, when I simulate the two cases, the generator angle theta differs by 180 degrees, so there are apparently two solutions to the problem.

3. if the system is linear, only one stationary solution exists; but the two systems shown here stay both in steady state (all constant values), but with different values from each other

See above, there are different solutions depending on theta.

I know that I used bad values for start values of the machine (UStartm UPhaseStart, PStart, QStart), chosen I don't remember how. I left them there on purpose, and chose not to make a load flow to find better initial values, to show that there is a problem.

If you provide the right power flow parameters, the library is designed to ensure that nonlinear solvers will converge to the corresponding solution. If you don't, your mileage may vary. There's nothing we can do about this - Achilles' heel of equation-based modelling when you have nonlinear implicit equations is that you may have convergence issues.

The only way to avoid this is to start from a reliable power flow and compute start values for all the nonlinear variables from that - this ensures that Newton's algorithm is initialized in the proximity of the steady-state power flow and converges to it, see the discussion in our ​paper. That's what we do in this library.

If you know any other way to guarantee the convergence to the "right" solution for a generic power grid initialization problem, that we can encode in Modelica in a declarative way, please let me know. I don't know any.

In fact, the points 1..3 are valid even if I give bad values for the initialization process: if only an initial point exists, the solver should either find it or not converge; In the unfortunate case it finds a bad initial point, the subsequent transient should not show all constant values.

No, because of theta, see above. It may be that one of the two solutions is stable and the other is unstable, but if you initialize at equilibrium with over 12 digits of precision, the system will stay close to the unstable equilibrium for quite a while. Unfortunately there is no way to tell the Newton solver to only converge to the stable solution, because it has no idea of stability. The ​HELM method can do this, but I'm not sure we can encode it in Modelica, and it is patented anyway :)

Well, maybe putting forward an issue in this way is not good for developers, which need simpler MWE's Maybe I can close the ticket for now, until when I find a simpler example to show the issue.

Before closing the ticket, I would suggest you to read the power flow values from the simulation results of Dymola and put them in the parameters. My bet is that then you will get the same results in both cases. If that is the case, you may claim that Dymola is a better tool than OpenModelica at finding the "right" solution when bogus power flow parameters are provided, and I may concede to that, though it may just be a matter of good vs. bad luck. Unfortunately, Dymola's methods are not published, contrary to OMC's, so I can't say why this happens.

@casella, just one more question: all the initialisation techniques in PowerGrids assume that all inter-connected generators share the system frequency, isn't it?

Yes. In fact, we further assume that the initial frequency is fixed to the nominal value, 50 or 60 Hz, as explained in the Tutorial

[massimo ceraolo]

Thank, you, Francesco, for your unpaid, thoughtful analysis.
I insisted not to follow the natural path (first load-flow, then transient), because I was convinced that one solution must exist.
Indeed when we power engineers connect a generator to a grid through a transformer, once the field voltage and input power are known, the operation point is unique.

My reasoning was partially wrong when I mentioned that just one solution must exist. That would have been true if the generator voltage were specified in amplitude and phase: i.e., it would be a phasor.
In the initialisation of my models (let us consider for a moment an isotropic machine, where Xd=Xq), indeed it is specified in terms of just amplitude, the phase being indirectly defined from the input power info. That makes the system different from an AC circuit, and therefore we are not sure we have one solution.

So, it should be not a big surprise finding that a second initial solution (with a different voltage from the first) exists. I expect, however, as you suggest, that this second solution is unstable, simply because, in real systems such as the one of this ticket, just one voltage can occur. I have no argument to state that a sufficiently stable numerical algorithm, when applied to a system, can or cannot converge to an unstable point, so I provisionally accept your idea that it can.

You also made a second logical step: that if a simulation starts from an unstable point, it stays stably on that point. That is harder to accept to me.
Thas would mean that while the equilibrium in the physical system is locally unstable when a numerical algorithm (say dassl) is applied to it, it is locally numerically stable, at least for small perturbations (consider that we can see from the plots some fluctuations of the solution, e.g. looking at GEN.port.VPu with OMEdit that currently shows every fluctuation, even tiny).

All this said, I resorted to the standard procedure to starting from a load flow solution.
You can see the results in Test2.mo => as expected the transient starts from initial machine voltage 1.0 and stays there the whole transient.

So we can conclude that:

• in nearly all PowerGrids systems it is very advisable to start from load flow analysis (with the possible exception of just a single generator feeding constant impedance load)
• otherwise it can happen that the starting point is on an unwanted initial point, possibly unstable in the physical system
• the fact that the found initial point stays stably constant during subsequent transient probably does not guarantee that the initial point is stable.

Now, I want to finish this ticket's saga with a different note.
I created this system to check on time-domain the behaviour of the PowerGrids SM under short circuit, with the transformer terms included as is common when time-domain simulations (with sine waves details) are made. It seems to me that no tests of that behaviour were included in the tests and examples of PowerGrids, probably because the focus is on phasor simulations and phasor simulations are normally run neglecting transformer terms.

So, if the models fom Test2 are run for 4s instead of 0.1s, one can see all Park's variables, plus three currents as a function of time.

I used the Entsoe machine, started at a no-load condition and nominal voltage, and then subjected to a three-phase short circuit. I checked the currents and torques as a function of time, comparing the results with those of a similar simulation made using a Matlab/Simscape model.
The results of this comparison are very good. Only very minor differences are seen.

[Francesco Casella]

Replying to ceraolo:

Thank, you, Francesco, for your unpaid, thoughtful analysis.

I'm trying to get prepared for the battle with domain-specific tools, it's good to have a tough sparring partner :)

I insisted not to follow the natural path (first load-flow, then transient), because I was convinced that one solution must exist.
Indeed when we power engineers connect a generator to a grid through a transformer, once the field voltage and input power are known, the operation point is unique.

The stable one, yes. Unfortunately there is no theoretical guarantee that you always converge to that.

My reasoning was partially wrong when I mentioned that just one solution must exist. That would have been true if the generator voltage were specified in amplitude and phase: i.e., it would be a phasor.
In the initialisation of my models (let us consider for a moment an isotropic machine, where Xd=Xq), indeed it is specified in terms of just amplitude, the phase being indirectly defined from the input power info. That makes the system different from an AC circuit, and therefore we are not sure we have one solution.

So, it should be not a big surprise finding that a second initial solution (with a different voltage from the first) exists. I expect, however, as you suggest, that this second solution is unstable, simply because, in real systems such as the one of this ticket, just one voltage can occur. I have no argument to state that a sufficiently stable numerical algorithm, when applied to a system, can or cannot converge to an unstable point, so I provisionally accept your idea that it can.

You can also try to linearize the model using OMC and then do an eigenvalue analysis to confirm that.

You also made a second logical step: that if a simulation starts from an unstable point, it stays stably on that point. That is harder to accept to me.

I never said that. I just said that if you start very, very close to it, it will take a very, very long time before the state trajectory eventually blows up.

Assuming that the unstable eigenvalue has the order of magnitude of a few radians per second, the doubling time of the positive exponential has the order of magnitude of a second. Hence, the distance from the equilibrium will double after about one second, and become 1000 times after 10 seconds (2¹⁰ = 1000), 1e6 times after 20 seconds, 1e9 times after 30 seconds, 1e12 times after 40 seconds. In fact, if you apply a 1e-4 amplitude perturbation to the excitation voltage at time = 1, the angle theta changes by a hundred degrees in 10 seconds.

All this said, I resorted to the standard procedure to starting from a load flow solution.
You can see the results in Test2.mo => as expected the transient starts from initial machine voltage 1.0 and stays there the whole transient.

So we can conclude that:

• in nearly all PowerGrids systems it is very advisable to start from load flow analysis (with the possible exception of just a single generator feeding constant impedance load)

Yes. This is very clearly stated in the tutorial, if you think it is not enough, we can strengthen the argument.

• otherwise it can happen that the starting point is on an unwanted initial point, possibly unstable in the physical system

Yes.

• the fact that the found initial point stays stably constant during subsequent transient probably does not guarantee that the initial point is stable.

Yes.

Now, I want to finish this ticket's saga with a different note.
I created this system to check on time-domain the behaviour of the PowerGrids SM under short circuit, with the transformer terms included as is common when time-domain simulations (with sine waves details) are made. It seems to me that no tests of that behaviour were included in the tests and examples of PowerGrids, probably because the focus is on phasor simulations and phasor simulations are normally run neglecting transformer terms.

So, if the models fom Test2 are run for 4s instead of 0.1s, one can see all Park's variables, plus three currents as a function of time.

I used the Entsoe machine, started at a no-load condition and nominal voltage, and then subjected to a three-phase short circuit. I checked the currents and torques as a function of time, comparing the results with those of a similar simulation made using a Matlab/Simscape model.
The results of this comparison are very good. Only very minor differences are seen.

That's very good to hear. I think we should include this test case in the Example package. We just need to categorize it properly. If you already made some other tests we could coordinate with Adrien, please let me know.

[massimo ceraolo]

Just in case it may be helpful I'll send privately to you a report on this check and the .slx file.