The backend does not compute strong components of initial equations correctly with if initial() equations

[Francesco Casella]

Please have a look at the PowerGrids.Examples.Tutorial.IslandOperation.TwoGeneratorsLocalInitialization model, using --tearingMethod=minimalTearing. I get the nonlinear initial system 472 containing 65 unknowns, belonging to GEN2, GRIDL2, LINE, and TGEN2. According to the way the model is designed, the initial equations should result in two separate strong components, one involving unknowns from GEN2, the other one involving GRIDL2, LINE, and TGEN2.

I manually went through all the equations using the debugger, and in fact I couldn't find any equation involving variables from both the subsets at the same time. So, unless I am missing something, I guess we not only have problems with nonlinear equation detection (see #5770), but also with Tarjan's algorithm and strong component detection.

Please note some of these equations are found inside if initial() clauses, this could be relevant to find the root cause of the issue.

@Karim, can you please check ASAP?

[Francesco Casella]

@Karim, I attached the library in a mail to you on 20/01/2020 22:12

[Andreas Heuermann]

It seems, that we are getting connect equations inside the Loop for GEN2 and TGEN2 which are not removed by removeSimpleEquations, maybe because they are with variables of a record (?)

44/44 (1): 0.0 = GEN2.GEN.port.v.im - TGEN2.vA.im   [unknown |0|0|0|0|]
55/55 (1): 0.0 = TGEN2.iA.im - GEN2.GEN.port.iGen.im   [unknown |0|0|0|0|]
57/57 (1): 0.0 = TGEN2.iA.re - GEN2.GEN.port.iGen.re   [unknown |0|0|0|0|]
58/58 (1): 0.0 = GEN2.GEN.port.v.re - TGEN2.vA.re   [unknown |0|0|0|0|]

This could be a clue on what's going on and why omc thinks it is one big loop.

[Francesco Casella]

@AnHeuermann, those equations show up in the regular set, and they are just fine. During initialization, they are replaced by other equations thanks to an if initial() clause, so those equations are not present in the initial equation set.

The problem is apparently due to the combination of initial equations and if initial() clause.

[Karim Adbdelhak]

If we have an equation like

  a = if initial() then x else y;

whe should just have

  a = y;

for the regular system and

  a = x;

for the initial one. Right?

Maybe we considering both inputs for causalization in both cases.

[Karim Adbdelhak]

There isn't any check for the initial operator, i think this is a potential fix. Needs some time to implement since the actual information if it belongs to an initial system is like four levels higher. Should take some time to change all necessary functions, i will do that tomorrow!

[Karim Adbdelhak]

I started a ​PR that implements rudimentary support for initial(). It did not seem to fix the problem. There could be done quite more eg:

• if-equation support (right now only if-expression)
• enhanced adjacency matrix support, information for tearing and index reduction. as far as i understood it goes wrong before those phases so i concentrated on fixing the regular adjacency matrix.
• convoluted usage of initial(). right now only plain call and negated variant. any convoluted version is treated as if it was a plain call.

I tried to find the actual calls of initial() in the loop equations to have more information but didn't find them. Francesco Could you provide an example equation? Do you think any of the points i mentioned could matter rn?

[Francesco Casella]

Replying to Karim.Abdelhak:

I started a ​PR that implements rudimentary support for initial().

So I understand it wasn't supported previously. I didn't know about that.

Using initial() in this way is indeed a bit tricky compared to the usual case, because normally you don't want regular equations to be ignored or replaced by others during initialization. However, this is absolutely essential in this case to replicate a traditional and robust initialization approach, which is used in domain-specific power system simulation codes, where of course one can write whatever he wants for initialization and simulation separately.

It did not seem to fix the problem. There could be done quite more eg:

• if-equation support (right now only if-expression)

In this case we really need if-equation support. However, for the time being you can restrict yourself to the case where the condition of the if-equation depends only on parameters having binding equations and fixe=true, and on the initial() operator. I don't think we'll ever need anything more than this.

• enhanced adjacency matrix support, information for tearing and index reduction. as far as i understood it goes wrong before those phases so i concentrated on fixing the regular adjacency matrix.

Getting the right incidence and adjacency matrices is essential

• convoluted usage of initial(). right now only plain call and negated variant. any convoluted version is treated as if it was a plain call.

What do you mean by "convoluted"? What I need here is to support logical expressions involving Boolean parameters and initial(). However, I guess the front end will take care of those parameters by evaluating them, since they are obviously structural parameters.

I tried to find the actual calls of initial() in the loop equations to have more information but didn't find them. Francesco Could you provide an example equation?

You find them in PowerGrids.Electrical.BaseClasses.OnePortAC, I report them here for your convenience:

equation
  if initial() and localInit == LocalInitializationOption.PV then
    // During local initialization, P,V is enforced at the connector towards
    // the grid, while Q,phase(V) is enforced at the port towards the component
    // Further two initial equations will be needed (e.g. in the controller)
    // to enforce P,V at the port or achieve actuator saturation if that is not possible
    CM.real(terminal.v*CM.conj(terminal.i)) = port.PStart;
    CM.'abs'(terminal.v) = port.VStart;
    port.Q = port.QStart;
    port.v.re*port.vStart.im = port.v.im*port.vStart.re;
  elseif initial() and localInit == LocalInitializationOption.PQ then
    // During local initialization, P,Q is enforced at the connector towards
    // the grid, while phase(V), phase(I) is enforced at the port towards the component
    // Further two initial equations will be needed (e.g. in the controller)
    // to enforce P,Q at the port or achieve actuator saturation if that is not possible
    terminal.v*CM.conj(terminal.i) = Complex(port.PStart, port.QStart);
    port.v.re*port.vStart.im = port.v.im*port.vStart.re;
    port.i.re*port.iStart.im = port.i.re*port.iStart.im;
  else
    // In all other cases, the port voltage and current are equal to
    // the connector voltage and current
    port.v = terminal.v;
    port.i = terminal.i;
  end if;

Based on the values of certain boolean parameters, certain equations are selected at initialization and possibly others during simulation.

The NF gets evaluates the parameters, so the flat model looks like this

  if initial() then
    Modelica.ComplexMath.real(PowerGrids.Types.ComplexVoltage.'*'.multiply(GEN1.GEN.terminal.v, Modelica.ComplexMath.conj(GEN1.GEN.terminal.i))) = GEN1.GEN.port.PStart;
    Modelica.ComplexMath.'abs'(GEN1.GEN.terminal.v) = GEN1.GEN.port.VStart;
    GEN1.GEN.port.Q = GEN1.GEN.port.QStart;
    GEN1.GEN.port.v.re * GEN1.GEN.port.vStart.im = GEN1.GEN.port.v.im * GEN1.GEN.port.vStart.re;
  else
    GEN1.GEN.port.v = GEN1.GEN.terminal.v;
    GEN1.GEN.port.i = GEN1.GEN.terminal.i;
  end if;

Note that the latter two equations involve Complex variables, so there are four scalar equations in each branch.

I guess handling this pattern is really trivial, just pick the former equations to build the initial and initial-lambda0 sets of equations, and pick the latter ones for the set of regular equations.

[Karim Adbdelhak]

Replying to Karim.Abdelhak:

I think i will try to provide some additional explanation:

• if-equation support (right now only if-expression)

The following is an if-equation;

if initial() then 
  x = a;
else
  x = b;
end if;

The same can be stated with an if-expression:

x = if initial() then a else b;

Right now we only support the second. There actually is a transformation phase from the first to the second since we need that for code generation anyways, but i am not quite sure when that happens. Note that the first allows multiple equations in each branch which can be beneficial as well as disadvantageous.

• enhanced adjacency matrix support, information for tearing and index reduction. as far as i understood it goes wrong before those phases so i concentrated on fixing the regular adjacency matrix.

Adjacency matrix: only plain occurrence for matching
Enhanced adjacency matrix: additional information about solvability
incidence matrix: Regularly used as a synonym for adjacency matrix but plainly wrong. we don't have any incidence matrices in the backend, only adjacency matrices.

e.g: z = y * f(x); with f(x) not invertable
Adjacency matrix entry: (z), (y),(x)
Enhanced adjacency matrix: (z, solvable), (y, nonlinear), (x, unsolvable)

• convoluted usage of initial(). right now only plain call and negated variant. any convoluted version is treated as if it was a plain call.

with convoluted i mean something like if b and initial(). we don't check for any stuff more than the plain call (and the negated version) in an if condition.

The PR got merged and i got another one with enhanced adjacency stuff pending so maybe we can see improvements in tonights tests.

[Karim Adbdelhak]

Replying to casella:

The NF gets evaluates the parameters, so the flat model looks like this

  if initial() then
    Modelica.ComplexMath.real(PowerGrids.Types.ComplexVoltage.'*'.multiply(GEN1.GEN.terminal.v, Modelica.ComplexMath.conj(GEN1.GEN.terminal.i))) = GEN1.GEN.port.PStart;
    Modelica.ComplexMath.'abs'(GEN1.GEN.terminal.v) = GEN1.GEN.port.VStart;
    GEN1.GEN.port.Q = GEN1.GEN.port.QStart;
    GEN1.GEN.port.v.re * GEN1.GEN.port.vStart.im = GEN1.GEN.port.v.im * GEN1.GEN.port.vStart.re;
  else
    GEN1.GEN.port.v = GEN1.GEN.terminal.v;
    GEN1.GEN.port.i = GEN1.GEN.terminal.i;
  end if;

Note that the latter two equations involve Complex variables, so there are four scalar equations in each branch.

I guess handling this pattern is really trivial, just pick the former equations to build the initial and initial-lambda0 sets of equations, and pick the latter ones for the set of regular equations.

Last night i was too tired, i missed that part! Okay as it seems we need to have initial() in if-equations handled and not only in if expressions! I added the enhanced adjacency matrix support and started a ​PR, but i will try to add this.

[Francesco Casella]

Replying to Karim.Abdelhak:

The following is an if-equation;

if initial() then 
  x = a;
else
  x = b;
end if;

The same can be stated with an if-expression:

x = if initial() then a else b;

In our case it will be more like

if initial() then
  a = b;
else 
  c = d;
end if;

which I guess becomes

0 = if initial() then (a-b) else (c-d)

In fact, we have two equations in each branch, as you see in comment:9, I'm not sure how this can be handled.

[Francesco Casella]

Replying to Karim.Abdelhak:

Last night i was too tired, i missed that part!

Sorry, I'm kind of a night owl.

Okay as it seems we need to have initial() in if-equations handled and not only in if expressions!

I think we basically only need to support this pattern

if initial() then
  // some equations
else
  // some other equations
end if;

everything else involving (structural!) parameters will be handled by the front-end.

Maybe there are also use cases where initial() is combined with time-varying conditions in conditional equations, but I'm not aware of them, I would not bother at this time.

I added the enhanced adjacency matrix support and started a ​PR, but i will try to add this.

OK, pls keep me posted.

[Karim Adbdelhak]

Replying to casella:

OK, pls keep me posted.

The old PR is merged, i made a ​new one with the if-equation changes. It does not seem to fix this particular problem, but maybe i am overlooking something.

Let's see what jenkins has to say about it and i will try to analyze what is happening.

[Karim Adbdelhak]

Replying to casella:

So I understand it wasn't supported previously. I didn't know about that.

Also a short comment on this one. For code generation everything was fine, we just lacked correct symbolic processing for matching/sorting/tearing/index reduction. Instead of only considering either the then or else case, depending on which system we analyze, we just took both.

Most of the time that would not cause crucial errors since some slightly too big loops don't usually cause a big loss in performance (sometimes it is even wanted e.g. DAEMode). Besides the operator itself does not seem to be used that often.

[Francesco Casella]

Replying to Karim.Abdelhak:

So I understand it wasn't supported previously. I didn't know about that.

Also a short comment on this one. For code generation everything was fine, we just lacked correct symbolic processing for matching/sorting/tearing/index reduction. Instead of only considering either the then or else case, depending on which system we analyze, we just took both.

Most of the time that would not cause crucial errors since some slightly too big loops don't usually cause a big loss in performance (sometimes it is even wanted e.g. DAEMode). Besides the operator itself does not seem to be used that often.

I understand. Unfortunately, in this case it is really crucial.

[Francesco Casella]

I just added a MWE that demonstrates the problem. The regular equations have a strong component with five unknowns, which is wrongly preserved during initialization, where in could be split into several smaller strong components.

[Karim Adbdelhak]

Thanks that minimal model helped a lot!

I realized that at some point the initial() call gets simplified to true for the initial system and false for the simulation system (after causalization). Seems great, but there is no simplification afterwards and the traverser just did not check if the condition of if expressions and equations can be simplified to true or false and just traversed both branches. I added that check in ​PR745 and it should work now.

Hopefully it speeds up some stuff and does not break anything!

[Francesco Casella]

The PR is merged, I just started a nightly build, I'll check that out myself ASAP.

Thanks!

[Francesco Casella]

The test case still doesn't work, and due to #5805 I can't understand why, I need to check from the command line.

However, I checked the structural information about the initial system, and it seems to me that the equations are now correctly split based on the if initial() condition.