Incorrect handling of algorithms with arrays by the backend

[casella]

@kabdelhak, while investigating #6199 and #6207, I stumbled upon something weird related to how algorithm sections involving arrays are handled by the backend.

I could reproduce the issue (or a related one) with a much simpler model in the context of #7815. Consider this test package:

package TestAlgorithmArray

  model EquationBased
    parameter Integer N = 2;
    parameter Real C = 1;
    parameter Real G = 1;
    Real T[N](each start = 1, each fixed = true);
    Real Ta[N], Tb[N];
    Real Qa[N], Qb[N];
  equation
    Ta[1] = 0;
    Tb[N] = 0;
    Qa[1] = 2*N*G*(Ta[1] - T[1]);
    Qb[N] = 2*N*G*(Tb[N] - T[N]);
    for i in 1:N-1 loop
      Qb[i] = N*G*(T[i+1] - T[i]);
      Qa[i+1] = -Qb[i];
      Tb[i] = (T[i+1] + T[i])/2;
      Ta[i+1] = Tb[i];
    end for;
  equation
    for i in 1:N loop
      der(T[i]) = N/C*(Qa[i] + Qb[i]);
    end for;
  end EquationBased;
  model AlgorithmBased
    parameter Integer N = 2;
    parameter Real C = 1;
    parameter Real G = 1;
    Real T[N](each start = 1, each fixed = true);
    Real Ta[N], Tb[N];
    Real Qa[N], Qb[N];
  algorithm
    Ta[1] := 0;
    Tb[N] := 0;
    Qa[1] := 2*N*G*(Ta[1] - T[1]);
    Qb[N] := 2*N*G*(Tb[N] - T[N]);
    for i in 1:N-1 loop
      Qb[i] := N*G*(T[i+1] - T[i]);
      Qa[i+1] := -Qb[i];
      Tb[i] := (T[i+1] + T[i])/2;
      Ta[i+1] := Tb[i];
    end for;
  equation
    for i in 1:N loop
      der(T[i]) = N/C*(Qa[i] + Qb[i]);
    end for;
  end AlgorithmBased;
end TestAlgorithmArray;

The two models should obviously be equivalent; they also give the same result in Dymola. However, for some reason, when I try to compile the AlgorithmBased model in OMC, I get

[1] 16:39:40 Translation Error
Internal error IndexReduction.pantelidesIndexReduction failed! System is structurally singular and cannot be handled because the number of unassigned equations is larger than the number of states. Use -d=bltdump to get more information.

[2] 16:39:40 Translation Error
Internal error Transformation Module PFPlusExt index Reduction Method Pantelides failed!

I suspect that a large number of the issues that are still plaguing the Buildings library are related to this. Can you please check?

[kabdelhak]

On first glance i see a major error. The adjacency rows for the algorithm seem to be faulty. I made an even smaller example:

package arr_eq_test
model eq
  Real T[2];
equation
  T[1] = 1;
  for i in 1:1 loop
    T[i+1] = T[i]*cos(time);
  end for;
end eq;

model arr
  Real T[2];
algorithm
  T[1] := 1;
  for i in 1:1 loop
    T[i+1] := T[i]*cos(time);
  end for;
end arr;

end arr_eq_test;

The adjacency matrix entry for the algorithm block in arr_eq_test.arr only shows entries for T[1] and not T[2], it seems like the iterator i+1 is not parsed correctly.

Seems like a crucial error and i don't know how that one slipped by for so long. I will see if i can fix it.

[casella]

Seems like a crucial error and i don't know how that one slipped by for so long. I will see if i can fix it.

Good that we caught it, eventually, and good that you have hopefully fixed it already. Let's see what happens with the testsuites.

[kabdelhak]

There is only one model failing, but while analyzing it i just realized that there could be something way more basic wrong here.

We talked about arrays in algorithms once and as far as i remember we came to following conclusion: If one element of an array is assigned in an algorithm, the whole array is assumed to be. Is that correct or does this only hold for initialization @casella ?

If that is the case it would be way easier to resolve these issues.

[casella]

There is only one model failing, but while analyzing it i just realized that there could be something way more basic wrong here.

We talked about arrays in algorithms once and as far as i remember we came to following conclusion: If one element of an array is assigned in an algorithm, the whole array is assumed to be. Is that correct or does this only hold for initialization @casella ?

It is correct. Modelica Specification Section 11.1.2:

Whenever an algorithm section is invoked, all variables appearing on the left hand side of the assignment operator := are initialized

• A continuous-time variable is initialized with its start value (i.e. the value of the start attribute).
• A discrete-time variable v is initialized with pre(v).
• If at least one element of an array appears on the left hand side of the assignment operator, then the complete array is initialized in this algorithm section.

This means that any given array must be completely handled in one algorithm section. You can't have two separate algorithm sections assigning values to different elements of an array. And if you only assign one (or a few) elements of the array explicitly, all the other ones are automatically assigned their start or pre() values.

I hope this help to make things simpler.

[kabdelhak]

• If at least one element of an array appears on the left hand side of the assignment operator, then the complete array is initialized in this algorithm section.

This only refers to initialization right? Or does this still hold for simulation?

[kabdelhak]

This particular problem is solved with #7853 but i will wait for the regression reports and see if it breaks anything.

Hopefully it also fixed something for the other issues.

[casella]

This only refers to initialization right? Or does this still hold for simulation?

I understand this always applies. There is no mention in the specification of the initialization phase. BTW, the reason why the non-discrete variables are always initialized to their start attribute also during simulation is precisely to avoid introducing some kind of "memory" that would not be appropriate, as explained in the spec.

[casella]

It passed the compiler testuite. Let's see what happens to the library testsuite.

[casella]

We did get some improvements in Buildings, see report, but not in the models reported in #6199 and 6207, I guess you need to dig more.

There were also some regressions in the test report, I think it's some frontend stuff, but I'm not sure.

[casella]

See comment here

[casella]

@kabdelhak you can use this MWE as a reference:

model M
  Real x[2](start = {3, 4}), y, z;
algorithm
  x[2] := y + z;
  x[2] := 2*x[2];
equation
  z = 2;
  x[2] = 10;
end M;

According to the Modelica Specification Section 11.1.2, I understand this is how this means.

• the algorithm assigns 3 to x[1] because of the start value rule
• the algorithm assigns x[2] to y + z, and then reassign 2x[2] to x[2]. This means that x[2] = 2(y+z), but the compiler shouldn't try to find an equivalent declarative formulation of the algorithm, because in general this may not even be theoretically possible (because of Turing's undecidability), so this is only between you and me :)
• so, the algorithm section is equivalent to

  equation
  x = f(y,z);
  function f
  input Real y; // RHS variable is a input
  input Real z; // RHS variable is an input
  output Real x[2]; // LHS variable is an ouput 
  algorithm
  // same algorithm as in the model + LHS initialization
  x := {2,3};
  x[2] := y + z;
  x[2] := 2*x[2];
  end f;

  Now comes matching. We have these equations:

  1.  z = 2;
  2.  x[2] = 10;
  3.  x = f(y,z);

  so obviously the first equation is matched with z, the second is matched with x[2], and the third is matched with whatever's left, i.e., x[1] and y. Note that in general algorithm sections are not matched with their LHS variables, as this example makes very clear.

Next comes causalization

z := 2;
x[2] := 10;
// implicit equation: unknowns are x[1] and y
res = x - f(x,y);

The solution is:

x[1] = 3;
x[2] = 10;
y = 3;
z = 2;

Note that y is solved "backwards", compare to the execution flow of the algorithm. No big deal, in this case, since the algorithm is linear, one Newton iteration would do, otherwise you'd need to iterate until convergence. There's no need to solve it symbolically.

I checked with Dymola, that runs M smoothly and produced the expected results. OMC 1.18.0 gives

[1] 14:31:05 Translation Error
Internal error IndexReduction.pantelidesIndexReduction failed! System is structurally singular and cannot be handled
because the number of unassigned equations is larger than the number of states. Use -d=bltdump to get more information.

[2] 14:31:05 Translation Error
Internal error Transformation Module PFPlusExt index Reduction Method Pantelides failed!

I'm not sure what happens with the latest nightly, you can check that.

[phannebohm]

What about the following model:

model M
  Real x[2](start = {3, 4}), y, z;
algorithm
  x[2] := y;
  x[2] := 5*z;
equation
  z = 2;
  x[2] = 10;
end M;

As I understand, y can take any value. Does that make the model singular? Do we need to detect that?

EDIT: The second assignment is already redundant, so this is definitely singular, right?

[kabdelhak]

The problem with the regressions was as follows. We have a tensor variable which gets assigned in a for loop, see small example:

model M
  parameter Integer m = 10;
  Real x[m,m,m];
algorithm
  for i in 1:m loop
    for j in 1:m loop
        for k in 1:m loop
            x[i,j,k] := i*sin(j*time+k);
        end for;
    end for;
  end for;
end M;

If we remove the subscripts from x it leads to a lot of element lookups. In general it tries hash table lookups for each element each time any element is assigned on the LHS. If we do not remove the subscripts it only looks up one element at a time, i think you can imagine how this scales so badly. In this case: WITHOUT FIX: ~24s WITH FIX: ~5s (scaling horribly for the unfixed version with rising m)

The fix is quite simple: Do not remove subscripts if they are constant. This seems weird because we do not have any constant subscripts here, but there are flattening steps for code generation which generate each variable and algorithm statement individually (if that is actually necessary is a whole other story, but that is how it currently works) and for those flattened versions this blows up terribly.

Bad byproduct: Models like the one presented where we use constant subscripts in the model itself still do not work because of the rule i just implemented. This is not really worrying for most cases, because generally this only occurs for variables which are iterated with a for loop, but conceptually this is a problem.

@phannebohm @casella any ideas on how to make it better without changing simcode/code generation entirely? For now it at least seems to work.

[casella]

@kabdelhak, I'm not aware of the details of how code generation works when algorithms are involved. I suspect there may be something fishy, but you need to explain to me how it works. Maybe we can set up a short meeting for that?

[casella]

@phannebohm

What about the following model:

model M
  Real x[2](start = {3, 4}), y, z;
algorithm
  x[2] := y;
  x[2] := 5*z;
equation
  z = 2;
  x[2] = 10;
end M;

As I understand, y can take any value.

That's correct.

Does that make the model singular?

Definitely.

Do we need to detect that?

I don't think so. According to my previous comment, this model will have to be solved numerically. A Jacobian will be computed, either numerically or symbolically by means of AD, and alas it will turn out to be singular at runtime. That's none of our business, if a modeller writes implicit algorithms in models and they turn out to be singular, that's his problem 😄

In other words, in this simple case you can tell by inspection that there is a structural singularity. However, in general algorithms can be arbitrarily complex and have unpredictable behaviour if there are iterations and conditional branches, so there's no point trying to catch these singularities at compile time. You may be able to do that in some special cases, but it's not worth spending time on this, we have a lot of higher priority work to do. The structural analysis assumes that x = f(y,z), if then function f() actually turns out to be singular at runtime, too bad.

EDIT: The second assignment is already redundant, so this is definitely singular, right?

I'm not really sure what you mean by "the second assignment is already redundant", but anyway the model is definitely singular. In fact, I tried it with Dymola and it simulated with y = 0, but IMHO this is questionable, one should at least get a warning that the model is not well posed.

[casella]

I'm not into the details of what the backend currently does with algorithms, but what I understand from these comments is that perhaps it may be trying to get too deep into the structural analysis of the algorithm itself and, as a consequence, may often end up down some endless rabbit hole (see the 130+ seconds spent by removeSimpleEquations, or this comment, or even fail, as #6199 and #6207 demonstrate.

I think the backend should really just see each algorithm section as a black box:

(v1, v2, ..., vn) = f_alg(w1, w2, ..., wn)

where the v's are the LHS variables (or arrays, which is also nice for the new backend, because we don't have to scalarize anything), the w's are the RHS variables other than the v's, and f_alg is a function with the w's as inputs, the v's as outputs, and the same algorithm of the original algorithm section, with an additional section at the beginning for LHS variable initialization to their start or pre values. That's it. One may add some optimization, e.g., by avoiding the initialization of variables that are assigned in the algorithm to the result of some expression that only contain already computed stuff, but that's not strictly necessary. In general, if one wants more sophisticated or deeper symbolic analysis, simplification and optimization, then he or she should use an equation section, not an algorithm section.

@perost, in fact, we may enact this view already in the frontend, by automatically generating those functions and replacing the algorithm sections by equations calling those function, as shown above. Then, we'd get rid of algorithm sections in models, and just need a backend capable of handling equations and functions. What do you think?

EDIT: if some variables among the v's are discrete, then we need to pass pre(v) to f_alg so they can be initialized to their pre() value. This is not necessary for continuous variables, because the frontend knows what the start value of v is an can substitute directly in the function algorithm body.

The equivalence is more complex in case of when statements in the algorithm, see later comment

[phannebohm]

Do we need to detect that?

I don't think so. According to my previous comment, this model will have to be solved numerically. A Jacobian will be computed, either numerically or symbolically by means of AD, and alas it will turn out to be singular at runtime. That's none of our business, if a modeller writes implicit algorithms in models and they turn out to be singular, that's his problem smile

In other words, in this simple case you can tell by inspection that there is a structural singularity. However, in general algorithms can be arbitrarily complex and have unpredictable behaviour if there are iterations and conditional branches, so there's no point trying to catch these singularities at compile time. You may be able to do that in some special cases, but it's not worth spending time on this, we have a lot of higher priority work to do. The structural analysis assumes that x = f(y,z), if then function f() actually turns out to be singular at runtime, too bad.

Good point!

EDIT: The second assignment is already redundant, so this is definitely singular, right?

I'm not really sure what you mean by "the second assignment is already redundant", but anyway the model is definitely singular. In fact, I tried it with Dymola and it simulated with y = 0, but IMHO this is questionable, one should at least get a warning that the model is not well posed.

Sorry I wasn't clear. I just mean the second assignment in the algorithm block x[2] := 5*z; together with the equations

equation
  z = 2;
  x[2] = 10;

If they are conflicting, this would be detected at runtime because the solver can't converge. But if the equations are in agreement, as they are in this example, any value of y is accepted without ever iterating, which I guess happens with Dymola. Of course this should produce a singular Jacobian in the first place, as you pointed out, and that should lead to a warning or even an error, since wild values of y may otherwise be used elsewhere.

[kabdelhak]

Looks like i was able to fix the FFT regressions: regression report Not quite sure if the speedups are also my merit.

Maybe we should still have a meeting about some general stuff @casella @phannebohm (@AnHeuermann ?). My points to discuss are:

• What kind of singularities do we need to detect?
• Invertability of Algorithms
• Behavior while solving algorithms in an algebraic Loop
• Can we treat constant accessed and assigned array variables as scalar variables (e.g. a[2] :=)? -> I suppose the answer is no due to the specification so instead: how do we deal with those without changing code gen entirely

[casella]

I'm not really sure what you mean by "the second assignment is already redundant", but anyway the model is definitely singular.

Sorry I wasn't clear. I just mean the second assignment in the algorithm block x[2] := 5*z; together with the equations

equation
  z = 2;
  x[2] = 10;

OK, now I get it, sorry. Singular square problems can have either infinitely many solutions (with x[2] := 5*x) or no solution at all (e.g. with x[2] := 4*x).

As a modeller, I think singular models are not well posed, even more so if they have infinitely many solutions, so they should be reported as flawed and abort the simulation. Apparently DS has other opinions.

If they are conflicting, this would be detected at runtime because the solver can't converge. But if the equations are in agreement, as they are in this example, any value of y is accepted without ever iterating, which I guess happens with Dymola.

Yes. I guess they are using some generalized solver using Moore-Penrose's pseudoinverse or something like that to get a solution in that case, and I guess it's the one with minimum norm (y = 0). While that is fine mathematically speaking, if the mathematical model describes something real, giving me only one solution when infinitely many are possible is highly questionable.

Of course this should produce a singular Jacobian in the first place, as you pointed out, and that should lead to a warning or even an error, since wild values of y may otherwise be used elsewhere.

I agree.

[casella]

Looks like i was able to fix the FFT regressions: regression report

Good!

Not quite sure if the speedups are also my merit.

They don't hurt, anyway.

Maybe we should still have a meeting about some general stuff @casella @phannebohm (@AnHeuermann ?). My points to discuss are:

* What kind of singularities do we need to detect?

IMHO, none, if algorithms are involved.

* Invertability of Algorithms

IMHO only numerically. If you can do it symbolically in some cases, that's fine, but don't waste too much time on that.

* Behavior while solving algorithms in an algebraic Loop

My proposal is to consider the algorithm section as a single equation with a function call. Unfortunately this doesn't cover algorithms with when clauses, I'm not sure this is within the scope of this ticket.

* Can we treat constant accessed and assigned array variables as scalar variables (e.g. `a[2] := `)?
  -> I suppose the answer is no due to the specification so instead: how do we deal with those without changing code gen entirely

You need to explain me how the current code gen works, then.

What about next Monday after the devmeeting?

[casella]

Some more thoughts about the case of when statements. Simple example case:

algorithm
  when v > 2 then
     x := y;
  end when;
  when v < 4 then
     x := -y;
  end when;

How do we turn this into an equation section? One option is to define a boolean variable (or array if there are multiple conditions) for each condition, bring them out in the equation section, and substitute them with if-statements in the algorithm, e.g.:

equation
  $cond1 = edge(v > 2); // Triggers event when v exceeds 2
  $cond2 = edge(v < 4); // Triggers an event when v goes below 4
  x = f_alg(y, pre(x), $cond1, $cond2);

function f_alg
  input Real y; // RHS variable in the algorithm
  input Real $pre_y; // pre() value of LHS variable
  input Boolean $cond1; // Triggering condition of first when-statement
  input Boolean $cond2; // Triggering condition of second when-statement
  output Real x; // LHS variable in the algorithm
algorithm
  x := $pre_x; // Initialization, v is a discrete variable because it is assigned in a when-section
  if edge($cond1) then
     x := y;
  end when;
  if edge($cond2) then
     x := -y;
  end when;

Similar handling applies for if-statements, of course

[casella]

Updated version of proposal to handle algorithm sections in the NF

• For each when-statement in an algorithm section, auxiliary boolean variable definitions (or arrays in case of multiple conditions) Boolean | Boolean[:] $cond_n = edge(<condition_n>); are added to the variable definitions section
• For each if-statement in an algorithm section, auxiliary boolean variable definitions Boolean $cond_n = <condition_n> are added to the variable definitions section
• Each algorithm section is substituted by an equation section with one function call to compute all the LHS variables
• The function has all RHS variables, pre() values of LHS variables, and event-triggering conditions as inputs, all LHS variables (or arrays) as outputs, and the same algorithm as the algorithm section, except that
  • Statements are added at the beginning of the algorithm to initialize all LHS variables to their start attributes (if continuous) or to their pre()-value (if discrete)
  • each when <condition_n> statement in the original algorithm section is replaced by if $cond_n in the function algorithm
  • each if <condition_n> statement in the original algorithm is replaced by if $cond_n in the function algorithm
• If it is (easily!) possible to determine that the initialization value of a certain variable will be always overwritten, the initialization statements may be omitted for optimization purposes.

  
  Boolean | Boolean[:] $cond1 = edge(<condition of first when statement>);
  Boolean | Boolean[:] $cond2 = edge(<condition of second when statement>);
  ...
  Boolean $cond_k = <condition of first if statement>;
  Boolean $cond_k+1 = <condition of second if statement>;
  ...
  equation
  (v1, v2, ..., vn) = f_alg_n(w1, w2, ..., wn, pre(vd1), pre(vd2), ..., $cond1, $cond2, ...)

function f_alg_n // one for each algorithm section input w1 "First RHS variable or array"; input w2 "Second RHS variable or array"; ... input pre_vd1 "pre-value of first discrete LHS variable"; input pre_vd2 "pre-value of second discrete LHS variable"; ... input Boolean | Boolean[:] $cond1 "First even-triggering condition"; input Boolean | Boolean[:] $cond2 "Second event-triggering condition"; ... output v1 "First LHS variable or array"; output v2 "Second LHS variable or array;" ... algorithm // Initialization section v1c := v1c.start; // First LHS continuous variable v2c := v2c.start; // Second LHS continuous variable ... v1d = pre_v1d; // First LHS discrete variable v2d = pre_v2d; // Second LHS discrete variable ... // The original algoritm section gets copied here with the following replacements: // when => if $cond_n // if => if $cond_n end f_alg_n;

I guess in case the original algorithm section contains variables with dot notation, e.g. a.b.c, the variable definitions in the function algorithm should be redefined as a_b_c or something like that.

[adrpo]

I guess you would also need to pass package constants as inputs. Also functions used in the algorithm section should be visible to the created function.

[casella]

I guess you would also need to pass package constants as inputs.

In fact, the function is defined locally in the model, so it "sees" the package constants exactly as the algorithm section does. The NF just needs to constant evaluate them before transforming the algorithm in a function.

Also functions used in the algorithm section should be visible to the created function.

Ditto.

@adrpo, in fact I understand the main consequence of your remark is that we shouldn't perform this transformation as the last step of NF or first step of BE, i.e., after flattening, but that in fact we should be perform it before flattening, hence in the NF. In this way, lookup of package constants and functions and whatever else is needed in the algorithm will works by construction, since the function is defined within the model that contains the algorithm section, so whatever is visible from the model is also visible from within the function. Once that has been done, standard flattening applies.

@perost, what do you think?

[phannebohm]

Some more thoughts about the case of when statements. Simple example case:

algorithm
  when v > 2 then
     x := y;
  end when;
  when v < 4 then
     x := -y;
  end when;

How do we turn this into an equation section? One option is to define a boolean variable (or array if there are multiple conditions) for each condition, bring them out in the equation section, and substitute them with if-statements in the algorithm, e.g.:

equation
  $cond1 = edge(v > 2); // Triggers event when v exceeds 2
  $cond2 = edge(v < 4); // Triggers an event when v goes below 4
  x = f_alg(y, pre(x), $cond1, $cond2);

function f_alg
  input Real y; // RHS variable in the algorithm
  input Real $pre_y; // pre() value of LHS variable
  input Boolean $cond1; // Triggering condition of first when-statement
  input Boolean $cond2; // Triggering condition of second when-statement
  output Real x; // LHS variable in the algorithm
algorithm
  x := $pre_x; // Initialization, v is a discrete variable because it is assigned in a when-section
  if edge($cond1) then

Why the double edge, is that on purpose? The defining equation for $cond1 already has that.

      x := y;
    end when;
    if edge($cond2) then
       x := -y;
    end when;

Similar handling applies for if-statements, of course

I like the approach. But one problem I see is if you have a when condition that depends on a LHS of the algorithm:

algorithm
  x := 2*y;
  when x > 1 then
    z : = 1;
  end when;

In this case it could be difficult to extract the condition. We could substitute the LHS with its RHS and continue with your approach, but that could lead to extremely long expressions. And I'm not sure if that's always possible. The part before the when statement might contain other conditional behavior.

[kabdelhak]

We have a bunch of general problems when converting algorithms to functions (Modelica Spec v3.6-dev 12.2 ):

A function cannot contain calls to the Modelica built-in operators der, initial, terminal, sample, pre, edge, change, reinit, delay, cardinality, inStream, actualStream, to the operators of the built-in package Connections, to the operators defined in chapter 16 and chapter 17, and is not allowed to contain when-statements.

I don't know if it is really beneficial to convert algorithms to functions because there are so many special cases to handle. On the new backend branch we added Input/Output lists to the algorithm record, maybe that suffices for efficient backend analysis.

[casella]

We have a bunch of general problems when converting algorithms to functions (Modelica Spec v3.6-dev 12.2 ):

A function cannot contain calls to the Modelica built-in operators der, initial, terminal, sample, pre, edge, change, reinit, delay, cardinality, inStream, actualStream, to the operators of the built-in package Connections, to the operators defined in chapter 16 and chapter 17, and is not allowed to contain when-statements.

And therein lies the rub 😢

I don't know if it is really beneficial to convert algorithms to functions because there are so many special cases to handle.

Of course not. I would then consider my proposal as a semi-formal description of what the backend should actually do, including handling all those things that are forbidden in functions.

What is essential is that we stick to the provisions of the language specification, in particular:

• An algorithm section is treated as an atomic vector-equation, which is sorted together with all other equations. For the sorting process (BLT), every algorithm section with N different left-hand side variables, is treated as an atomic N-dimensional vector-equation containing all variables appearing in the algorithm section. This guarantees that all N equations end up in an algebraic loop and the statements of the algorithm section remain together (Section 11.1.2).
• An algorithm section specifies how to compute output variables as a function of given input variables. A Modelica tool may actually invert an algorithm section, i.e., compute inputs from given outputs, e.g., by search (generate and test), or by deriving an inverse algorithm symbolically (Section 11.1.

@kabdelhak are we arleady doing this? The problem posted at the beginning of this long thread now works, and that fixed some models in the library testsuite, but this MWE

model M
  Real x[2](start = {3, 4}), y, z;
algorithm
  x[2] := y;
  x[2] := 5*z;
equation
  z = 2;
  x[2] = 10;
end M;

still fails with the latest nightly build, so I understand the current backend does not fully comply with the specification.

On the new backend branch we added Input/Output lists to the algorithm record, maybe that suffices for efficient backend analysis.

OK. However, we should try hard to see if there is some reasonable fix that we can apply to the existing backend to get Buildings going. It may not be efficient (that would be for the new backend) but it should at least work.

[casella]

More specifically, I checked again the three test cases that are reported as exemplary of the 52 Buildings 7.0.1 models failing in the backend:

• Buildings.Fluid.HeatExchangers.ActiveBeams.Examples.CoolingAndHeating, ticket #6207
• Buildings.Experimental.DistrictHeatingCooling.Validation.HeatingCoolingHotWaterSmall, ticket #6200
• Buildings.Experimental.DistrictHeatingCooling.Plants.Validation.Plant_Carnot_T_ClosedLoop, ticket #6199

All three cases have structural issues related to either causalization or initialization. The flattened models reveals that in all three cases there are some initial algorithms, and also some algorithms in the regular section in the first case. I guess that is no coincidence 😄

I'm not sure if the reason why they fail is the same as in the MWE of the previous comment, i.e. that you are forcing the matching to the LHS variables of algorithms section, and that does not allow to complete the causalization/index reduction step properly, or if there is any other reason for that. However, I think it's worth exploring a fix that would make the MWE work, and see if those issues get resolved.

Let's discuss this next Monday after the devmeeting.

[casella]

Further analysis of #6207 reveals that we have the following code:

algorithm
  beaCooHea.conHea.mod.temDif_mod.i := 1;
  for j in 1:2 loop
    if beaCooHea.conHea.mod.temDif_mod.u > beaCooHea.conHea.mod.temDif_mod.xd[j] then
      beaCooHea.conHea.mod.temDif_mod.i := j;
    end if;
  end for;
  beaCooHea.conHea.mod.temDif_mod.y := Buildings.Utilities.Math.Functions.cubicHermiteLinearExtrapolation(beaCooHea.conHea.mod.temDif_mod.u, beaCooHea.conHea.mod.temDif_mod.xd[beaCooHea.conHea.mod.temDif_mod.i], beaCooHea.conHea.mod.temDif_mod.xd[beaCooHea.conHea.mod.temDif_mod.i + 1], beaCooHea.conHea.mod.temDif_mod.yd[beaCooHea.conHea.mod.temDif_mod.i], beaCooHea.conHea.mod.temDif_mod.yd[beaCooHea.conHea.mod.temDif_mod.i + 1], beaCooHea.conHea.mod.temDif_mod.dMonotone[beaCooHea.conHea.mod.temDif_mod.i], beaCooHea.conHea.mod.temDif_mod.dMonotone[beaCooHea.conHea.mod.temDif_mod.i + 1]);

so the algorithm first computes and index and then uses it to access some data in order to compute some other output. You need to scalarize variable subscripts on the RHS

[casella]

The other two test cases only have initial algorithms, but they only contain asserts, so the root cause must be different.

[casella]

Further analysis of the other two cases reveals that in both ones it involves redundant initial equations on the air side of heat exchangers, which features volumes that are directly connected to the atmosphere (fixed pressure sinks). This is particularly clear in model HeatingCoolingHotWaterSmall, where if -d=bltdump is turned on, a singularity involving only two equations is spotted, where the two equations are p = p_source and p = p_start.

So, it seems that #6200 is correlated to #6199, but debugging is probably much easier in #6200.

It may be that the redundant initial equation is not handled because one of the two has a parameter that may be changed at runtime. We need to check what happens if we set Evaluate = true or if we put some modifiers with literals on that parameter, maybe that solves the problem, in which case we should make the algorithm less picky and declare that parameter as structural in this case.

EDIT: the following test case works in OMC 1.19.0-dev

model M2
  Real x(stateSelect = StateSelect.prefer), y(stateSelect = StateSelect.prefer), z, w1;
  Real X(stateSelect = StateSelect.avoid),Y(stateSelect = StateSelect.avoid);
  parameter Real x_start = 1;
  parameter Real y_start = 1;
  parameter Real C = 2;
  parameter Real c = 1;
equation
  X = C*x;
  Y = C*y;
  der(X) = z + w1;
  der(Y) = -z;
  x = y;
  x = c;
initial equation
  x = x_start;
  y = y_start;
end M2;

so it seems the problem is different

[casella]

Plan: @kabdelhak tries to fix #6207, @casella tries to see if setting p_start to constant is enough for the redundant initial equation algorithm to work.