Simple equation involving Complex numbers is solved incorrectly

[Francesco Casella]

I have a power system model containing the following equation

vAt = k*vA;

where vAt and vA are SI.ComplexVoltage variables and k is a Complex constant. In this specific case, k is a Real constant quantity, so

k = Complex(r, 0.0);

where r is a Real parameter.

In the specific model I am trying to run, vAt is computed first, so OMC tries to solve this equation for vA. This should be indeed trivial to do symbolically, and should lead to this result:

vA.re := DIVISION(vAt.re,r);
vA.im := DIVISION(vAt.im,r);

What I get instead is

$cse5 := Complex.'constructor'.fromReal(r, 0.0);
linear torn system with 2 unknowns, tearing variable is vA.re
(torn) vA.im := DIVISION($cse5.re * vA.re - vAt.re, $cse5.im)
(residual) $cse5.im * vA.re + $cse5.re * vA.im - vAt.im = 0

Unfortunately OMC does not solve this equation symbolically, which should be trivial in this case, as k.im is known to be a literal zero constant.

Even worse, it tries to solve this as an implicit system of linear equations, and it picks a tearing variable which leads to a division by k.im, which is in fact a literal zero, so it causes a singularity:

(torn) vA.im = (r*vA.re - vAt.re)/0.0 = 0.0/0.0

I don't think we should try to fix the tearing here, but rather try to make sure that equations dealing with Complex numbers are really brought down to their scalar components, so that all symbolic simplifications can be applied, particularly in those cases where the equations are actually decoupled, as in this case.

I tried to reproduce the problem with this simple model

model foobar
  Complex vAt;
  Complex vA;
  parameter Complex k = Complex(r,0.0);
  parameter Real r = 10;
equation
  vAt = k*vA;
  vAt = Complex(time,2*time);
end foobar;

but of course by Murphy's law, this time OMC selects vA.im as the tearing variable, so it can actually solve the equation. In any case it does so in this way;

(torn) vA.re := DIVISION(vAt.re + k.im * vA.im, k.re)
(residual) k.im * vA.re + k.re * vA.im - vAt.im = 0

which is not the right way to go for the reasons discussed above.

Anyone can help me to make sure that the foobar model is solved as

vA.re := DIVISION(vAt.re,r);
vA.im := DIVISION(vAt.im,r);

I guess that would also solve my original problem, because no tearing will be used then.

Thanks!

[Vitalij Ruge]

I guess a smart rule for solve complex equation in Backend would be difficult (relat to #4354).
But, it's possible to symbolically solve all n x n linear systems with +maxSizeSolveLinearSystem=m where n <=m . So +maxSizeSolveLinearSystem=2 work well for foobar.

[Francesco Casella]

I tried in my actual model and still got an error:

assert | debug | <p>division by zero at time 0, (a=2001.028986583865) / (b=-0),
where divisor b expression is: __OMC__2$QR$h$1

I traced back how __OMC__2$QR$h$1 is computed and apparently it should be non-zero, but unfortunately I cannot check the values of already computed variables and parameters when debugging.

[Francesco Casella]

In the meantime a workaround is to set --maxSizeLinearTearing=0, but we need to solve this issue for good.

[Vitalij Ruge]

I had the train of thought:
for the simplication:

vA.re := DIVISION(vAt.re,r);

I would prefer

final parameter Complex k = Complex(r,0.0);

so k.im should replace 0.

[Vitalij Ruge]

forbar with final final parameter Complex k = Complex(r,0.0);:

+d=dumpindxdae:

Equations (2, 2) 
======================================== 
1/1 (1): vA.im = vAt.im / k.re   [dynamic] 
2/2 (1): vAt.im = 2.0 * time   [dynamic] 

Equations (2, 2) 
======================================== 
1/1 (1): vA.re = vAt.re / k.re   [dynamic] 
2/2 (1): vAt.re = time   [dynamic] 

[Francesco Casella]

Replying to vitalij:

I had the train of thought:
for the simplication:

vA.re := DIVISION(vAt.re,r);

I would prefer

final parameter Complex k = Complex(r,0.0);

Brilliant! For the test case I posted, you are absolutely right, in order to be able to perform the symbolic simplification, final is actually mandatory.

However, this still doesn't solve the issue with the original model, where the binding for k is actually found in the equation section

  parameter SI.PerUnit rFixed "Fixed transformer ratio";
equation
  vAt = k*vA;
  k = Complex(rFixed);

[Vitalij Ruge]

by remove prefix parameter for k:

model foobar
  Complex vAt;
  Complex vA;
  Complex k = Complex(r);
  parameter Real r = 10;

equation
  vAt = k*vA;
  vAt = Complex(time,2*time);

end foobar;

I get the SimCode dump:

        b-vector: 
        7: vA.re=DIVISION(vAt.im - k.re * vA.im, k.im) [Real ] 
        8: k.re * vA.re + (-k.im) * vA.im - vAt.re (RESIDUAL) 

with devision by zero k.im = 0.

[Vitalij Ruge]

Complex k = Complex(r,0);

work well !!

removeSmpleEquation can work here!
Complex k = Complex(r) != Complex(r, 0)
because Complex(r) is a function call. Possible inliningComplex.'constructor'.fromReal would be a great idea.

[Vitalij Ruge]

model from comment:8 start to run if inline for Complex.'constructor'.fromReal will be fixed

e.g. by add a function body ( no empty algorithm section, otherwise inlining is failed) for Complex.'constructor'.fromReal in Complex3.2.2.mo

      output Complex result "Complex number";
    algorithm
      result := Complex(re=re, im=im) "Complex number";

[Francesco Casella]

@vitalij, thanks for this further analysis. In fact, I assumed that Complex(r) was just a shorthand notation for Complex(r,0) and that they would lead to the same behaviour (of course they should...).

I'm not sure why the MSL code of this constructor was built using binding equations on the output rather than with a one-liner algorithm, but I guess that was to make inlining even simpler. I guess we should handle that as well, changing the MSL code as you suggest would rather be a hack.

So, it seems like the root cause of the problem is the lack of inlining of Complex constructor with one argument and the other one defaulting to zero. Can someone take care of that?

[Vitalij Ruge]

I guess I can fix inline at the weekend. (For empty function-body)

[Vitalij Ruge]

Note:
I guess tearing selection is relate to #3511, because OM don't use the FunctionTree and maybe Complex.'constructor'.fromReal arguments don't lead a dependence of k.im.

[Vitalij Ruge]

maybe fixed with ​PR2310

[Francesco Casella]

Replying to vitalij:

maybe fixed with ​PR2310

Unfortunately as of v1.13.0-dev-573-gf10312d57, the foobar model is still solved by means of 2x2 linear equations and tearing.

[Vitalij Ruge]

I guess the model need some of the following modeficatio:

  // ----------- 1 ------------
    final parameter Complex k = Complex(r,0.0);
  // ----------- 2 ------------
     Complex k = Complex(r,0.0);
  // ----------- 3 ------------
      protected
      parameter Complex k = Complex(r,0.0);

otherwise the modelica tool can't (because not allowed) set k.im =0 and remove the loop.

[Vitalij Ruge]

@casella Anyway, I guess you original model should work, now. Can you try it, please?

[Francesco Casella]

Sorry, of course I need to add the final qualifier or remove the parameter declaration and make it a variable with binding.

model foobar1
  Complex vAt;
  Complex vA;
  final parameter Complex k = Complex(r,0.0);
  parameter Real r = 10;
equation
  vAt = k*vA;
  vAt = Complex(time,2*time);
end foobar1;

model foobar2
  Complex vAt;
  Complex vA;
  Complex k = Complex(r,0.0);
  parameter Real r = 10;
equation
  vAt = k*vA;
  vAt = Complex(time,2*time);
end foobar2;

Both work as expected now, ending in an all-assignment model.

My original test model also works.

Thanks @vitalij!

[Francesco Casella]

Unfortunately this test case still doesn't work properly

model foobar3
  Complex vAt;
  Complex vA;
  Complex k = Complex(r);
  parameter Real r = 10;
equation
  vAt = k*vA;
  vAt = Complex(time,2*time);
end foobar3;

Apparently, Complex(r) is still treated differently from Complex(r,0.0). Of course the workaround is to add the second input to the constructor, but you can't always do that, if you are using somebody else's library (e.g. the MSL).

@vitalij, can you give one more look to this test case also? Thanks!

[Martin Sjölund]

The easy fix would be to fix the inlining so it works even if the algorithm section is empty :)

[Vitalij Ruge]

fixed with ​PR2331
improved test with ​PR909

[Francesco Casella]

Thanks @vitalij, everything runs fine now