Post optimization for lambda0 system failing for large power grid model

[AnHeuermann]

New ticket from https://trac.openmodelica.org/OpenModelica/ticket/5805#comment:16.

@casella is working on a bigger model that uses homotopy to initialize, but the backend is failing at calling BackendDAEUtil.postOptimizeDAE in ​Initialization.mo#L262, so the resulting lambda0 system is empty.

The example is not public yet, so only available via mail.

[casella]

Probably also related to #5806

[AnHeuermann]

I was on the wrong track to begin with.
Apparently we have two different functions to generate our SimCode object. One for ODE-mode (createSimCode) and one for DAE-mode (generateModelCodeDAE).

While createSimCode has initialEquations_lambda0 the dae one will always be empty. And for our example we are using daeMode, but I wasn't aware of that.
That was the cause for a lot of confusion for me.

And with that a model that is using daeMode will do nothing for lambda=0. That's not so good for our examples.

---

So do we need an initial system for lambda=0 for DAE-mode?
Probably, but let's try without first.

First of all let's set
int (*functionInitialEquations_lambda0)(DATA *data, threadData_t*); from our C callback to be NULL if initialEquations_lambda0 in SimCode.SIMCODE is empty.
The let our C runtime handle such a case by using the normal initialization system with lambda=0, e.g. calling functionInitialEquations instead of functionInitialEquations_lambda0.

[AnHeuermann]

I added some workaround for this to ​https://github.com/OpenModelica/OpenModelica/pull/7095.

The title of this ticket is very misleading. The problem I suspected to be the cause for #5805 is simply a totally different one.

The next step should be to add the initialEquations_lambda0 system to SimCode for dae models as discussed in #5806.

[casella]

Replying to AnHeuermann:

While createSimCode has initialEquations_lambda0 the dae one will always be empty. And for our example we are using daeMode, but I wasn't aware of that.
That was the cause for a lot of confusion for me.

Sorry about that, I should have been more explicit. These grid models need a unique combination of factors to work really well:

• daeMode simulation, as they are sparse as such, but dense if causalized into odeMode
• use of sparse solver for initialization, otherwise it simply takes forever for models of moderate or large size
• use of homotopy for initialization, to reduce the number of nonlinear variables that need a start value, and particularly to avoid them entirely in controllers, which are fully linear when simplified
• no tearing of initial equations, to avoid throwing away precious start values of nonlinear variables and instead use meaningless start values of linear variables that may eventually influence nonlinear variables in the torn section. For a better solution, see #6108

Regarding lambda0, I'm still unsure that it is necessary for power grid models. I know from experience that detailed power plant models require it, because the structural simplifications that are made possible by the simplified equations are essential to get convergence of the initial equations when lambda = 0, and the initialization solver fails abjectly without them.

From this point of view, notice that even though the equations with structural simplifications are equivalent to those without, but with numerical lambda = 0, meaning that they have the same solution, they may not be equally easy to solve by Newton's equation. This only became clear to me recently. For example, sqrt(x) = 10 and x = 100 are obviously equivalent if x.min = 0, but while the former is strongly nonlinear, so it can suffer a lot from a bad start value, the latter is linear, so it doesn't bother at all about it. So, if I had to solve

homotopy(sqrt(x) + 0.01*x, sqrt(x)) = 10

dealing with lambda = 0 separately allows to reduce the equation to

x = 100

which is trivially solved. Not doing it leads to

sqrt(x) + lambda*0.01*x = 10

which is hareder to process symbolically and is significantly nonlinear even when lambda = 0 if kept as such.

So do we need an initial system for lambda=0 for DAE-mode?
Probably, but let's try without first.

Sounds good that we try. Just bear in mind that we have evidence that handling lambda = 0 separately and exploiting the corresponding structural simplifications can do wonders in other cases.

First of all let's set
int (*functionInitialEquations_lambda0)(DATA *data, threadData_t*); from our C callback to be NULL if initialEquations_lambda0 in SimCode.SIMCODE is empty.
The let our C runtime handle such a case by using the normal initialization system with lambda=0, e.g. calling functionInitialEquations instead of functionInitialEquations_lambda0.

This sounds good. However, is it really different from what is happening now? I understood that the initial-lambda0 case actually only set the start values to the unknowns, and then the homotopy started at lambda = 0, which should be exactly the same as you're trying to do.

Do I miss something?

[AnHeuermann]

From this point of view, notice that even though the equations with structural simplifications are equivalent to those without, but with numerical lambda = 0, meaning that they have the same solution, they may not be equally easy to solve by Newton's equation.

I agree. But I would like to close this issue and continue on implementing this in #5806.

---

This sounds good. However, is it really different from what is happening now? I understood that the initial-lambda0 case actually only set the start values to the unknowns, and then the homotopy started at lambda = 0, which should be exactly the same as you're trying to do.
Do I miss something?

Before ​https://github.com/OpenModelica/OpenModelica/pull/7095 was merged for lambda=0 nothing was happening. So start values are set and then the homotopy process started with lambda=0.333 (or so). That is of course even harder then to start with the actual*lambda+simplified*(1-lambda) for lambda=0.

[casella]

Replying to AnHeuermann:

From this point of view, notice that even though the equations with structural simplifications are equivalent to those without, but with numerical lambda = 0, meaning that they have the same solution, they may not be equally easy to solve by Newton's equation.

I agree. But I would like to close this issue and continue on implementing this in #5806.

Sounds good.

---

This sounds good. However, is it really different from what is happening now? I understood that the initial-lambda0 case actually only set the start values to the unknowns, and then the homotopy started at lambda = 0, which should be exactly the same as you're trying to do.
Do I miss something?

Before ​https://github.com/OpenModelica/OpenModelica/pull/7095 was merged for lambda=0 nothing was happening. So start values are set and then the homotopy process started with lambda=0.333 (or so). That is of course even harder then to start with the actual*lambda+simplified*(1-lambda) for lambda=0.

OK, I see. I'll try out the nightly build on my grid tests and report if it works better.

Thanks!

[casella]

Just tried it with Type1_N_4_M_4 test case. Newton solver hangs forever when trying to solve the large nonlinear system during the lambda = 0 iteration. Besides, no LOG_NLS_V output for that system, alas, so it is impossible to figure out what's going on. You can try it yourself with the privately sent test package.

[casella]

Problem solved, see ticket:5806#comment:12.

The fact that the lambda = 0 step failed (without symbolic simplifications) maybe was due to #6348, we can continue the investigation using @christiankral's test case there.