Missing inlining of initial equations causes inefficiencies with Complex initial equations

[Francesco Casella]

Please consider the following test case (also attached)

loadString("
model M1
  parameter Real P = 1000;
  parameter Real Q = 2000;
  parameter Real V_re = 100;
  parameter Real V_im = 100;
  parameter Real I_re(fixed = false);
  parameter Real I_im(fixed = false);
initial equation
  P = V_re*I_re + V_im*I_im;
  Q = V_re*I_im - V_im*I_re;
end M1;

model M2
  import CM = Modelica.ComplexMath;
  parameter Real P = 1000;
  parameter Real Q = 2000;
  parameter Real V_re = 100;
  parameter Real V_im = 100;
  parameter Complex I(re(fixed=false), im(fixed=false));
initial equation
  Complex(P,Q) = Complex(V_re, V_im) * CM.conj(I); 
end M2;
");getErrorString();

setCommandLineOptions("-d=dumpSimCode");
loadModel(Modelica);
buildModel(M1);getErrorString();
buildModel(M2);getErrorString();

M2 is completely equivalent to M1, only written down in a more compact and convenient way, by exploiting Complex records and function overloading. There should be no performance penalty when M2 is simulated instead of M1, otherwise one would prefer writing M1 instead of M2, which is against the whole concept of high-level declarative modelling.

This is the output from the back-end for M1

5:  (LINEAR) index:0 jacobian: true
	variables:
index: -1: I_im (no alias)  initial: 	no arrCref index:() [] 
	b-vector:
	1: I_re=DIVISION(P - V_im * I_im, V_re) [Real ]
	2: V_re * I_im + (-V_im) * I_re - Q (RESIDUAL)
	Jacobian idx: 0
	3: I_re.$pDERLSJac0.dummyVarLSJac0=DIVISION((-V_im) * I_imSeedLSJac0, V_re) [Real ]
	4: $res.1.$pDERLSJac0.dummyVarLSJac0=V_re * I_imSeedLSJac0 - V_im * I_re.$pDERLSJac0.dummyVarLSJac0 [Real ]

this is the output of the back-end for M2:

1:   $TMP_Complex1 := Complex(V_re * I.re + V_im * I.im, V_im * I.re - V_re * I.im);
	
2: P=$TMP_Complex1.re [Real ]
3: Q=$TMP_Complex1.im [Real ]

With M1, OMC recognizes the intial equations as a linear system, and solves it via tearing. In fact, the way this is done is a bit questionable, because the solution would become unnecessarily singular when V_re = 0, which is perfectly legitimate, but that's another story.

With M2, OMC does not recognize this fact, and eventually a nonlinear iterative solver would be used. Of course I understand that in this specific case this is no big deal, as the Newton algorithm will converge in just one step.

My point is that with the current handling of Complex equations, in more complicated cases OMC may miss opportunities for useful symbolic handling, which could instead be critical.

I would recommend that Complex equations are symbolically transformed into Real equations, so that all kind of optimizations can be performed in the same way as with the flat formulation of the problem with Real variables. In this specific case, M1 and M2 would eventually be handled exactly in the same way by the solver.

One way to do so could be that in the special case of equations involving records with only two Real components r1 and r2

<lhs> = <rhs>

they are transformed into a pair of Real equations

<lhs.r1> = <rhs.r1>
<lhs.r2> = <rhs.r2>

instead of introducing auxiliary variables of any kind.

I reckon this can be trivially implemented in MetaModelica.

[anonymous]

Maybe the description should read as follows?

<lhs.r1> = <rhs.r1>
<lhs.r2> = <rhs.r2>

[Francesco Casella]

Replying to anonymous:

Maybe the description should read as follows?

<lhs.r1> = <rhs.r1>
<lhs.r2> = <rhs.r2>

Absolutely :)

[Francesco Casella]

Milestone changed to 1.13.0 since 1.12.0 was released

[Francesco Casella]

The problem with this model is the handling of initial equations. In fact, for regular equations

model M3
  import CM = Modelica.ComplexMath;
  Real P = 1000;
  Real Q = 2000;
  Real V_re = 100 + sin(time);
  Real V_im = 100 + cos(time);
  Complex I;
equation
  Complex(P, Q) = Complex(V_re, V_im) * CM.conj(I);
end M3;

the back-end eventually figures out that the equation is a linear one and solves it as such.

I guess the reason why the reported problem happens for initial equations but not for regular ones is that ​initOptModules does not include function inlining, which is required to further break down the result of the complex multiplication in the rhs into its real and imaginary components.

@lochel, @wbraun, @vwaurich, is there any good reason why function inlining is not included in the initOptModules? Could you please add that and make it default?

[Willi Braun]

I do not see any good reason to avoid inlining of function for the initial equation system.
I can try to add it and we will see, if the testsuite maybe has something against.

[Willi Braun]

Ah, okay. Actually for the correct handling of the regular equations in M3 it's not inlineFunction responsible, and not even the Backend. It is already done by the Frontend, but the Frontend does handle the initial equations different for some reasons.

loadString("
model M3
  import CM = Modelica.ComplexMath;
  Real P = 1000;
  Real Q = 2000;

  Real P_p = 1000;
  Real Q_p = 2000;

  Real V_re = 100 + sin(time);
  Real V_im = 100 + cos(time);

  parameter Real V_re_p = 100;
  parameter Real V_im_p = 100;

  Complex I;
  parameter Complex I_p(re(fixed=false), im(fixed=false));
initial equation
  Complex(P_p,Q_p) = Complex(V_re_p, V_im_p) * CM.conj(I_p);
equation
  Complex(P, Q) = Complex(V_re, V_im) * CM.conj(I);
end M3;
");getErrorString();

setCommandLineOptions("-d=dumpdaelow,dumpindxdae,dumpinitialsystem,dumpSimCode");
loadModel(Modelica);
buildModel(M3);getErrorString();

=>

########################################
dumpdaelow
########################################
[...]
Equations (8, 8)
========================================
1/1 (1): P = 1000.0   [binding]
2/2 (1): Q = 2000.0   [binding]
3/3 (1): P_p = 1000.0   [binding]
4/4 (1): Q_p = 2000.0   [binding]
5/5 (1): V_re = 100.0 + sin(time)   [binding]
6/6 (1): V_im = 100.0 + cos(time)   [binding]
7/7 (1): Q = V_im * I.re - V_re * I.im   [dynamic]
8/8 (1): P = V_re * I.re + V_im * I.im   [dynamic]

[...]

Initial Equations (1, 2)
========================================
1/1 (2): Complex(P_p, Q_p) = Complex(V_re_p * I_p.re + V_im_p * I_p.im, V_im_p * I_p.re - V_re_p * I_p.im)   [dynamic]

Also technically Complex is a record and not a function, so inlining can't handle this currently, so we probably we need something else in the Backend for such purpose.
@adrpo or @perost
On the other side how will the new Frontend act here? Should this done by the Backend or Frontend?

However, it seems our testsuite doesn't have anything against adding lateInlineFunction to the initOptModules.
​https://test.openmodelica.org/hudson/view/PullRequestTesting/job/OpenModelica_TEST_PULL_REQUEST/5378/
@lochel Do you have something?
BTW: I've noted that wrapFunctionCall is also not used by the initialization. Is there any good reason as well?

[Francesco Casella]

I just re-ran the M3 test with v1.13.0-dev-1096-g6770e6300. The output of flattening M3 with the NF is

class M3
  Real P = 1000.0;
  Real Q = 2000.0;
  Real P_p = 1000.0;
  Real Q_p = 2000.0;
  Real V_re = 100.0 + sin(time);
  Real V_im = 100.0 + cos(time);
  parameter Real V_re_p = 100.0;
  parameter Real V_im_p = 100.0;
  Real I.re "Real part of complex number";
  Real I.im "Imaginary part of complex number";
  parameter Real I_p.re(fixed = false) "Real part of complex number";
  parameter Real I_p.im(fixed = false) "Imaginary part of complex number";
initial equation
  Complex.'constructor'.fromReal(P_p, Q_p) = Complex.'*'.multiply(Complex.'constructor'.fromReal(V_re_p, V_im_p), Modelica.ComplexMath.conj(I_p));
equation
  Complex.'constructor'.fromReal(P, Q) = Complex.'*'.multiply(Complex.'constructor'.fromReal(V_re, V_im), Modelica.ComplexMath.conj(I));
end M3;

so it seems to me that the NF handles the equation and the initial equation exactly the same way, unless something which is not shown in the flattened model (e.g., annotations) actually happens to be different.

If I use -d=optdaedump, the output of pre-optimization module normalInlineFunction (simulation), which is the first module to be run, is

1/1 (1): P = 1000.0   [binding |0|0|0|0|]
2/2 (1): Q = 2000.0   [binding |0|0|0|0|]
3/3 (1): P_p = 1000.0   [binding |0|0|0|0|]
4/4 (1): Q_p = 2000.0   [binding |0|0|0|0|]
5/5 (1): V_re = 100.0 + sin(time)   [binding |0|0|0|0|]
6/6 (1): V_im = 100.0 + cos(time)   [binding |0|0|0|0|]
7/7 (1): Q = V_im * I.re - V_re * I.im   [dynamic |0|0|0|0|]
8/8 (1): P = V_re * I.re + V_im * I.im   [dynamic |0|0|0|0|]

so it seems to me that inlining is actually applied also to the Complex constructor, which happens to contain annotation(Inline=true).

Conversely, the initial equations are never subject to an inlining pass, so the initial complex equation

1/1 (2): Complex(P_p, Q_p) = Complex(V_re_p * I_p.re + V_im_p * I_p.im, V_im_p * I_p.re - V_re_p * I_p.im)   [dynamic |0|0|0|0|]

treads through the backend unscathed.

I guess the fact the inlining (as well as wrapFunctionCall, as noticed by @wbraun) was missing on the initial equations hasn't been a problem so far, because the equations that needed inlining or CSE were always found in the simulation equations, e.g. for Media or MultiBody models, while initial equations are simple things such as x = x_start or der(x) = 0.

However, in the case of models using Complex numbers, it is very well possible that initial equations requiring this processing are used, as demonstrated by the M3 test case.

@wbraun, @lochel, can you please coordinate to add calls to normalInlineFunction and wrapFunctionCalls to the initial equations also? I guess this shouldn't take much effort.

Thanks!

[Lennart Ochel]

The pre-optimization modules are applied on the shared DAE used for simulation and initialization. And the pre-optimization module normalInlineFunction is apparently also traversing the initial equations, which means they should actually be inlined. Unfortunately, I have no time for an extensive investigation right now.

[Per Östlund]

Replying to casella:

Conversely, the initial equations are never subject to an inlining pass, so the initial complex equation

1/1 (2): Complex(P_p, Q_p) = Complex(V_re_p * I_p.re + V_im_p * I_p.im, V_im_p * I_p.re - V_re_p * I_p.im)   [dynamic |0|0|0|0|]

treads through the backend unscathed.

The equation is actually inlined already in BackendDAECreate.lowerEqn (ironically with a comment saying TODO: remove inline), i.e. in the very first phase of the backend where the backend DAE structure is created. When the normalInlineFunction module is run there are no function calls left to inline in the model.

What you see above is thus not calls to the Complex constructor, but Complex record expressions that unfortunately looks the same as function calls in the output. So it seems that the actual issue is that some part of the backend doesn't handle record expressions as well as it should, and not an inlining issue.

M3 fails to compile with the new frontend for some reason though, even though the optdaedump looks pretty much identical. One difference in the flat model is that the new frontend uses the correct constructor, while the old frontend incorrectly uses the default one (which happens to give the same result in this case). But they're all inlined anyway, so it shouldn't matter.

[Francesco Casella]

Replying to perost:
I'm not sure if I can follow all the details. My point is that, if you look at M3, there are an initial equation and a regular equation that look exactly the same

initial equation
  Complex(P_p,Q_p) = Complex(V_re_p, V_im_p) * CM.conj(I_p);
equation
  Complex(P, Q) = Complex(V_re, V_im) * CM.conj(I);

The only difference is that V_re_p, V_im_p, and I_p ar parameters instead of variables. Does that justify the different handling by the back-end, i.e., that the complex equation is not split into its scalar counterparts?

[Francesco Casella]

In fact, the basic issue I haven't yet understood is: what is the function in the backend that splits complex equations involving Complex record expressions into their scalar counterparts?

[Per Östlund]

I did some more investigation on why M3 fails with the new frontend, and it turns out my initial analysis was wrong regarding inlining when using the old frontend.

Both equations are completely inlined when using the new frontend, because it uses the fromReal constructor from Complex that has the Inline = true annotation. It actually gets inlined even if the annotation is removed, because BackendDAECreate.lowerEqn forces inlining for complex equations. But then it fails because createNonlinearResidualEquationsComplex is called, which has no case for Record = Record (it expects one side to be a function call when the other is a record expression).

The old frontend incorrectly selects the default constructor for Complex instead (which shouldn't be used when there's a userdefined constructor), which for some reason is marked as Inline = false. This means that they're never inlined by lowerEqn or any of the pre-optimization phases.

Fixing the inline flag for default constructors is trivial, just need to change the appropriate Inline.NO_INLINE to e.g. Inline.DEFAULT_INLINE in Static.elabCallArgs2. But this still won't make the constructors be inlined using the old frontend, because the Inline module has no handling for default record constructors (which uses a different data structure compared to normal functions). That should also be relatively easy to fix, since all that's needed is to take the function arguments and create a record expression from them.

But even if these two issues are fixed for the old frontend it will just run into the same issue as the new frontend with createNonlinearResidualEquationsComplex.

So to summarize: I can fix the old frontend so that the constructor calls gets inlined, but someone needs to fix the backend first so that it can handle the inlined equations correctly (which means fixing the backend so that it can handle M3 with -d=newInst).

[Francesco Casella]

@perost, based on your inputs here, I just opened #5236, which only involves plain record constructors, and reproduces the issue. I hope that makes it easier to fix the backend.

I don't think it's worth spending time fixing the old backend, since it has other issues with Complex numbers. Getting them to work right with the new front-end is fine for me.

As a side note, I was perhaps incorrectly labeling as "inlining" also the action of splitting record equations into scalar equations. It is true, as @lochel pointed out in yesterday's discussion, that inlining is carried out by the backend also on initial equations, since the right-hand-side constructor holds the result of the complex multiplication in the optdaedump.

However, for some reason the record = record equation is split into its scalar components by the back-end, while the record = record initial equation is not. This is the subject of #5236

[Francesco Casella]

Rescheduled to 1.14.0 after 1.13.0 releasee

[Francesco Casella]

After #5236 has been fixed, the output of the NF on M3 is

######################################## 
dumpdaelow 
######################################## 
 
Equations (8, 8) 
======================================== 
1/1 (1): Q = V_im * I.re - V_re * I.im   [dynamic |0|0|0|0|] 
2/2 (1): P = V_re * I.re + V_im * I.im   [dynamic |0|0|0|0|] 
3/3 (1): P = 1000.0   [binding |0|0|0|0|] 
4/4 (1): Q = 2000.0   [binding |0|0|0|0|] 
5/5 (1): P_p = 1000.0   [binding |0|0|0|0|] 
6/6 (1): Q_p = 2000.0   [binding |0|0|0|0|] 
7/7 (1): V_re = 100.0 + sin(time)   [binding |0|0|0|0|] 
8/8 (1): V_im = 100.0 + cos(time)   [binding |0|0|0|0|] 
 
Initial Equations (2, 2) 
======================================== 
1/1 (1): Q_p = V_im_p * I_p.re - V_re_p * I_p.im   [dynamic |0|0|0|0|] 
2/2 (1): P_p = V_re_p * I_p.re + V_im_p * I_p.im   [dynamic |0|0|0|0|] 

Therefore, when using the NF, the problem is solved. This will be the default behaviour in 2.0.0, and can already be achieved by setting -d=newInst