Processing arrays of records fails

[ceraolo]

This is similar but different from #4606
I hope the solution is as straightforward as it was in that case!
Consider the following code:

package NumEqBug
  model TestLineZb
    parameter Real y[:] = {8, 10} ;
    Complex Z[2, 2];
  equation
      Z =LineZb(
        y=y );
  end TestLineZb;
  function LineZb
    import Modelica.ComplexMath;
    input Real y[2];
    output Complex Z[2, 2];
  algorithm
    Z:=fill(Complex(1,1),2,2);
  end LineZb;
end NumEqBug;

Checking TestLineZb I get 4 equations and 8 variables. I would expect 8 equations and 8 variables.
Note that y[] was included in the code even if it is not used because otherwise in the flattened code the function call is not included.
Tested with
v1.13.0-dev-188-g2c5818d (64-bit) for Win

[perost]

Fixed in 193c439.

[perost]

Actually, this just fixes checkModel, but it seems like the backend is also wrong since it says 16 equations/8 variables when you try to compile the model. So I guess we should keep the ticket open until that's also fixed.

[casella]

Rescheduled to 1.14.0 after 1.13.0 releasee

[ceraolo]

The purpose here is to build a two-dimensional array of complex numbers.
A working solution for the same purpose is in comment:17 of ticket #4055, which uses loops instead of fill() function.

[casella]

If you make y a constant

package NumEqBug
    model TestLineZb
        constant Real y[:] = {8, 10} ;
        Complex Z[2, 2];
    equation
            Z =LineZb(
                y=y );
    end TestLineZb;
  
  function LineZb
    import Modelica.ComplexMath;
    input Real y[2];
    output Complex Z[2, 2];
  algorithm
    Z:=fill(Complex(1,1),2,2);
  end LineZb;
end NumEqBug;

the test case works. Apparently, there is a problem with functions taking arrays of scalars as inputs and record arrays as outputs. No idea why this should be the case.

[casella]

Karim, are you interested at having a look at this?

[casella]

Marking this as a blocker for 1.14.0, Complex numbers should work in that version.

[Karim.Abdelhak]

Since this still does not work i reopen it.

Somehow the dimension of the complex equation
Z = {{Complex(1.0, 1.0), Complex(1.0, 1.0)}, {Complex(1.0, 1.0), Complex(1.0, 1.0)}};
Which is generated by the NF, ends up being (4,4) instead of (2,2,2) . Probably a conversion problem from the NF structure to BE. I will be looking into this.

[adrpo]

@perost: is this an NF issue?

[perost]

Replying to adrpo:

@perost: is this an NF issue?

The types look correct when flattening with -t, i.e. both Z and the rhs has type Complex[2, 2]. But with records and functions there are types hiding everywhere, so there might be some type issues that shows up if e.g. the backend inlines the call.

Karim can investigate the issue, and if it turns out that the NF is setting the wrong type on something I can help fix it if needed.

[Karim.Abdelhak]

I found where the dimension gets computed wrong and fixed it, but now there is an issue in the code generation in SimCodeUtil:

[/home/kab/workspace/OpenModelica/OMCompiler/Compiler/SimCode/SimCodeUtil.mo:6458:7-6458:48:writable] 
Error: Internal error solving array equation: 1 : Z = NumEqBug.LineZb(y) for variable: Z[2,1].re

I looked into it and it fails because (Expression.isMatrix(e1) or Expression.isArray(e1))
where e1 is Z, fails. The actual variables Z[i,j].re and Z[i,j].im are not generated, i guess this is because it is a matrix of records (?). I don't know if this is due to wrong typing from the frontend or a backend issue.

Since pull requests are not possible right now i just changed line 1859 from the function lowerArrayEqn in BackendDAECreate.mo from
ds = List.map1(ds, intMul, i);
to
ds = listReverse(i::listReverse(ds));

Can someone have a look if all types regarding that are correct in the frontend?

[ceraolo]

I understand that this is the only remaining ticket on complex numbers among those listed in
#4835.

I just want to say that this issue comes from a very real-life application.

For me its solution is mandatory (or nearly so) if you want to use simulate multi-conductor electric lines in a quasi-stationary way.
This in turn is commonly done when simulating electrified railway lines, in which the timespan requires the quasi stationary approach, and lines are multi-conductor. I've already published several papers on this.

Thank everybody for dedicating time to this issue.

[perost]

Replying to Karim.Abdelhak:

Can someone have a look if all types regarding that are correct in the frontend?

I had a quick look at this, and the only place in SimCodeUtil where Expression.isMatrix(e1) or Expression.isArray(e1) is used is in createSingleArrayEqnCode.

The type of the expression doesn't matter in this case, what those functions check is what kind of expression e1 is. I.e. isArray will only be true if e1 is an Expression.ARRAY() (isMatrix will never be true for the NF, it doesn't create matrix expressions).

The types looks correct in this case though, both sides of the equation in Complex[2, 2] as expected.

[Karim.Abdelhak]

This is related to ticket #5350

We decided that functions which create arrays of records will not be implemented for the next release since it would take too much time. Fundamental changes for creating arrays of records have to be made in SimCode. This will be addressed later in the projects regarding better array processing.

For now we would have following suggestion:
Try to remove functions which have arrays of records as outputs as much as possible.
(prevented with -d=-nfScalarize maybe)

In this case the function does not get removed since the input is not constant, even if the input is not used. But the function call could be replaced by the fill operator (even if it depends on the input). The fill operator would afterwards also be scalarized, like it already would be.
Maybe detect functions which only have a single (array record) output and replace them in the model?

If this is easy to do, we should prevent functions of this kind (for now).

@perost @adrpo is this viable?

[ceraolo]

@Karim.Abdelhak, I'm not sure that I've understood well your post.
Do you propose for this case to remove the function call and put the fill in the calling modelica code instead of the function?
If this is what you propose I would recommend to not spend you time on this.

This ticket comes from my attempt of making all the models I've developed over the years in Dymola able to run also inside OM. So it is a long-term goal, and I don't have pressure for doing this fast. As far as I can use OM for my modelica models, I can mention this on my papers. That's the main reason for this ticket.

My models using functions creating complex matrices don't have just a fill() function inside, but a lot of computations instead: the function I reported in this ticket is just a minimal case of a much more complex real-life situation.

So, if I understand well, to have this issue solved I must wait for 2.0.0, a thing that I easily can do.

[Karim.Abdelhak]

Replying to ceraolo:

My models using functions creating complex matrices don't have just a fill() function inside, but a lot of computations instead: the function I reported in this ticket is just a minimal case of a much more complex real-life situation.

Okay in that case we should really generally work out a better concept on arrays of records in general.

So, I understand well, to have this issue solved, I must wait for 2.0.0, a thing that I easily can do.

As i mentioned already we will address this in future array research, i will postpone it for now!

[Karim.Abdelhak]

The computation for the size has been updated in ​PR313 but there are still open issues with general computation of arrays of records.

We figured out a sensible way of doing it and formulated a draft. It will change quite a lot in the backend, code generation and runtime, therefore it will take some time to implement it and actually get it running properly. We will keep this ticket open and post any updates on this issue.

[casella]

@Karim.Abdelhak, any news on this issue?

[casella]

Update, I ran this test case with a recent nightly build. Frontend works fine, backend works until end of preoptimization

Variables (8) 
======================================== 
1: Z[2,2].im:VARIABLE()  "Imaginary part of complex number"NumEqBug.TestLineZb, Complex type: Real [2,2] 
2: Z[2,2].re:VARIABLE()  "Real part of complex number"NumEqBug.TestLineZb, Complex type: Real [2,2] 
3: Z[2,1].im:VARIABLE()  "Imaginary part of complex number"NumEqBug.TestLineZb, Complex type: Real [2,2] 
4: Z[2,1].re:VARIABLE()  "Real part of complex number"NumEqBug.TestLineZb, Complex type: Real [2,2] 
5: Z[1,2].im:VARIABLE()  "Imaginary part of complex number"NumEqBug.TestLineZb, Complex type: Real [2,2] 
6: Z[1,2].re:VARIABLE()  "Real part of complex number"NumEqBug.TestLineZb, Complex type: Real [2,2] 
7: Z[1,1].im:VARIABLE()  "Imaginary part of complex number"NumEqBug.TestLineZb, Complex type: Real [2,2] 
8: Z[1,1].re:VARIABLE()  "Real part of complex number"NumEqBug.TestLineZb, Complex type: Real [2,2] 
 
 
Equations (1, 8) 
======================================== 
1/1 (8): Z = NumEqBug.LineZb(y)   [dynamic |0|0|0|0|] 

Then it fails on Pantelides

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

@Karim, can you estimate how much work is still left to get this to run?

[casella]

@ceraolo, we have an easy fix for this based on inlining. Would this be ok for your actual problem?

[Karim.Abdelhak]

Replying to casella:

@ceraolo, we have an easy fix for this based on inlining. Would this be ok for your actual problem?

The basic question is: Do you use functions with arrays of records that cannot be reduced to one line? In that case my fix would not work.

[ceraolo]

Replying to Karim.Abdelhak:

Replying to casella:

@ceraolo, we have an easy fix for this based on inlining. Would this be ok for your actual problem?

The basic question is: Do you use functions with arrays of records that cannot be reduced to one line? In that case my fix would not work.

This ticket is rather old for my research. So for my personal usefulness I could lower priority to high or even normal.

But for the general usefulness of OM it is still very important IMO. This ticket came from the need I had to automatically compute the electrical parameters of power transmission lines at t=0.
In case this is possible (and with my modelica programs is already possible with Dymola) the user can change, at each simulation, the geometry of the line (number of conductors, height, horizontal position etc.): the Modelica tool at t=0, using geometry data, computes electrical parameters, then starts simulation using the computed parameters.

This is a procedure computer programs specifically designed for power systems (e.g. EMTP, ATP) have been doing for decades. Nowadays modelica tools can do the same, but not (yet) OM.

Three years ago we at Pisa published a paper on this usage of Modelica. We published it taking as reference Dymola, but maybe there is still room for a second paper made based on OM, if this ticket is solved.

Unfortunately, computing complex matrices of lines requires a rather lengthy algorithm: it cannot be compressed into a single line

To summarise: I don't need this ticket to be solved soon. However, if it is solved I could try to write a paper exploiting the features attainable after this solution, that would add visibility to OM.

See also comment 17.

A final note. To attain the result described in this comment I had also to file a ticket to Dassault Systèmes, since had found three issues with functions outputting matrixes made of complex numbers in Dymola, that after a year were solved.

[casella]

@ceraolo, could you post one such function for us to see? I am wondering if you really need functions, or if you could do the same with an initial algorithm.

[ceraolo]

I enclose here one example. In case of need, I can privately share the whole package with you, along with the paper that describes its usage and usefulness.
I don't know whether it can be converted into an initial algorithm. But I'm sure that this conversion, if possible, would reduce dramatically the readability of the program.

function LineZ "Compute matrix of longitudinal impedance per metre of a multi-conductor line"
  import Modelica.Constants.*;
  import Modelica.ComplexMath;
  import Modelica.Utilities.Streams;
  input Integer n "number of conductors in line";
  input Real x[n] "horizontal abscissas of conductors (m)";
  input Real y[n] "vertical abscissas of conductors (m)";
  input Real r[n] "equivalent tubular radii (GMR's) (m)";
  input Real R1[n] "conductors lineic resistance (ohm/m)";
  input Real rho = 100 "ground resistivity";
  input Real f = 50 "frequency";
  output Complex Z[n, n] "Computed matrix (ohm/m)";
  output Complex Zself, Zmutual "Transposed line self and mutual impedance";
protected
  constant Complex j = Complex(0, 1) "Imaginary unit";
  Real D "Distance between a conductor and the image of another";
  Real d "Distance between two conductors";
  Real Theta "Theta angle for Carson's formulas";
  //  Real k_s=0.7788;
  //exp(-mu_r/4)=exp(-0.25)
  Real L "generic inductance";
  Real P "P coefficient from Carson's original paper";
  Real Q "Q coefficient from Carson's original paper";
  parameter Real a0 = 4 * pi * sqrt(5) * 1e-4 * sqrt(f / rho) "multiplies D in a-parameter from EMTPs Theory Book";
  //  a0*D=a; a is the same as "r" in Carson's paper and "csi" in Uribe's ipst03 paper
  // where D=2h for self impedanze, D=D'ik for mutual impedance
  Real R_g = pi ^ 2 * f * 1e-7;
  Real w = 2 * pi * f;
  Real a;
  //=a0*D oppure 2*a0*H; parametro adimensionale.
algorithm
//self impedances with ground return:
  for i in 1:n loop
    a := a0 * 2 * y[i];
    P := pi / 8 - sqrt(2) / 6 * a + a ^ 2 / 16 * (0.672784 + log(2 / a));
    Q := (-0.03860784) + 0.5 * log(2 / a) + sqrt(2) / 6 * a;
    L := mue_0 / (2 * pi) * log(2 * y[i] / r[i]) "L when soil is pure conductive";
    Z[i, i] := R1[i] + 4 * w * P / 1.e7 + j * w * (L + 4 * Q / 1.e7);
  end for;
//Carson (25) and (30)
//mutual impedances with ground return:
  for i in 1:n loop
    for jj in 1:i - 1 loop
      Theta := atan(abs((x[i] - x[jj]) / (y[i] + y[jj])));
      d := sqrt((x[i] - x[jj]) ^ 2 + (y[i] - y[jj]) ^ 2);
      D := sqrt((x[i] - x[jj]) ^ 2 + (y[i] + y[jj]) ^ 2);
      a := a0 * D;
      P := pi / 8 - sqrt(2) / 6 * a * cos(Theta) + a ^ 2 / 16 * (0.672784 + log(2 / a)) * cos(2 * Theta) + a ^ 2 / 16 * Theta * sin(2 * Theta);
      Q := (-0.03860784) + 0.5 * log(2 / a) + sqrt(2) / 6 * a * cos(Theta);
      L := mue_0 / (2 * pi) * log(D / d) "L when soil is pure conductive";
      Z[i, jj] := 4 * w * P / 1.e7 + j * w * (L + 4 * Q / 1.e7);
    end for;
  end for;
//distanza fra conduttore i e jj
//distanza fra conduttore i e immagine di jj
//nella seguente formula -0.0386+0.5*log(2/a) = 0.5(0.6159315-log(a)) del theory book
// in quando 0.5*0.6159315-0.5*log(2)=-0.03860784
//Carson (26) and (31)
  for i in 1:n loop
    for jj in i + 1:n loop
      Z[i, jj] := Z[jj, i];
    end for;
  end for;
  if n == 3 then
    Zself := 1 / 3 * (Z[1, 1] + Z[2, 2] + Z[3, 3]);
    Zmutual := 1 / 3 * (Z[1, 2] + Z[2, 3] + Z[3, 1]);
  end if;
  annotation(
    Documentation(info = "<html>
<p>BundleZ usa le formule di Carson, nell'approssimazione valida per a&LT;0.25.</p>
<p>In sostanza vva benissimo allla frequenza industraiale e d&agrave; risultati sostanzialmente identici a Line Constants.</p>
<p>Anche il BundleZ va abbastanza bene: la scelta di Paris era ben giustificata dalla disponibilit&agrave; di calcolo dell'epoca e 
dal fatto che lui operava con linee a notevole distanza dal suolo.</p>
<p>L&apos;uso della maggior precisione offerta qui &egrave; giustificata da:</p>
<p>- facilità; di calcolo odierna</p>
<p>- congruenza con LineConstants</p>
<p>- presenza nel caso dei sistemi ferroviari di conduttori al suolo (rotaie) o addirittura interrate (dispersore lineare di terra)</p>
<p>- migliore confrontabilit&agrave; con i risultati di Lamedica che usava per l&apos;appunto lineConstants.</p>
<p>Un confronto fra i risultati di BundleZb, BundleZ e Line Constants &egrave; riportato nel documento Carson Validation.docx.</p>
</html>"));
end LineZ;

[casella]

Of course we should support this kind of things eventually.

As a temporary workaround, one may put this very same code inside an initial algorithm.

@Karim, would that work, or would we have the same limitations at the end of the day?