Is static tearing of parameter-dependent linear systems a good idea?

[casella]

Maybe not. Consider this test case

loadString("
model M
  parameter Real a = 1;
  parameter Real b = 0.5;
  parameter Real c = 0.5;
  parameter Real d = 1;
  Real x,y;
equation
  a*x + b*y = 1;
  c*x + d*y = -1;
end M;
"); getErrorString();

setCommandLineOptions("-d=dumpSimCode");getErrorString();
simulate(M, simflags="-lv=LOG_STATS");getErrorString();
simulate(M, simflags="-lv=LOG_STATS -override=b=0");getErrorString();
simulate(M, simflags="-lv=LOG_STATS -override=b=1e-9");getErrorString();

The back-end solves the linear system by tearing like this:

6: y=DIVISION(1.0 - a * x, b) [Real ]
7: 1.0 + d * y + c * x (RESIDUAL)

which is fine when b=0.5, but is definitely not a good idea when b=0, and also not very numerically sound when b=1e-9.

My conclusion would be that we should by default avoid solving for variables that are multiplied by a parameter in the torn variable section.

If one wants the efficiency of full tearing, then he or she should turn on -d=evaluateAllParameters, which would at least avoid problem when b = 0 because that dependency will be totally eliminated.

Unfortuntely, that won't help much when b = 1e-9, unless the tearing heuristics somehow tries to avoid solving w.r.t. variables which are multiplied by a small coefficient.

Comments?

[ptaeuber]

The tearing heuristic would never choose the equation 6: y=DIVISION(1.0 - a * x, b) [Real ] to solve for y if b=0 because then it is treated as an impossible assignment. If you want to set b to zero you have to recompile the model and the other equation will be chosen as inner equation.

Besides that I changed the function which checks for parameters equal to zero from isZero to isZeroOrAlmostZero some time ago. Nevertheless, the threshold is chosen a bit too small (1e-15). I think changing it to 1e-6 should be fine?

It would definitively lead to much larger systems if we do not allow to divide by parameters.

[casella]

Replying to ptaeuber:

The tearing heuristic would never choose the equation 6: y=DIVISION(1.0 - a * x, b) [Real ] to solve for y if b=0 because then it is treated as an impossible assignment. If you want to set b to zero you have to recompile the model and the other equation will be chosen as inner equation.

I understand this. However, I have one test case from power systems where this would really not be an option (and also sound quite weird to the end user), see #4354. In that case, I'd get a singularity if V_re = 0, which is perfectly legitimate, as it implies a voltage phasor with an angle of +/- 90 deg. For that application, we have a specific requirement to be able to change the parameters at runtime, to avoid the huge code building times when running different simulations.

Besides that I changed the function which checks for parameters equal to zero from isZero to isZeroOrAlmostZero some time ago. Nevertheless, the threshold is chosen a bit too small (1e-15). I think changing it to 1e-6 should be fine?

You need to have a threshold that still guarantees reasonable rounding errors. 1e-6 would be ok in an ideal world where all variables and parameters are scaled to unity. But this is not necessarily the case. For instance, if the parameter in question is a power in a a large electrical power system, it may have order of magnitude 1e9, so 1e-6 would mean relative precision of 1e-15. We need a way to correctly take into account the scaling.

At the very least, the threshold should be 1e-6*p.nominal. I wonder if there could be some other criteria in case the nominal attribute is not set.

In any case, my example demonstrates, this doesn't really work if you change the parameter at runtime after compilation, which is always possible.

It would definitively lead to much larger systems if we do not allow to divide by parameters.

Of course. One recommmendation could be that if you get such a division by zero (or division by a small number, that would require a different macro in the code to catch the exception), you should get a meaningful error (for division by zero) or warning (for division by a small number) that points out that the root cause is linear tearing suggests a cure, e.g., to set some flag that avoids division by a parameterm, or at least to turn off tearing of linear systems.

With the current implementation, a "normal" user would just get a division by zero or a very numerically dirty solution without any clue so as to the cause and remedies. In the case of a numerically dirty solution, he may not even be aware of that, which is also bad.

[ptaeuber]

Ok, I see your point. Since there is already a flag -d=allowImpossibleAssignments, which allows to divide by almost everything, I suggest to add a config flag such as --tearingStrictness instead, which can take the values casual, strict and veryStrict, where veryStrict won't even allow dividing by parameters. I would leave the default to be strict.

For the normal strict tearing I will set the threshold to 1e-6*p.nominal for parameters and think about some kind of appropriate error message/warning in case something bad happens after overriding.

Maybe it would be enough to provide a warning that something bad could happen when a parameter is overridden with a small value or zero and tearingStrictness has not been veryStrict?

[casella]

Replying to ptaeuber:

Ok, I see your point. Since there is already a flag -d=allowImpossibleAssignments, which allows to divide by almost everything, I suggest to add a config flag such as --tearingStrictness instead, which can take the values casual, strict and veryStrict, where veryStrict won't even allow dividing by parameters. I would leave the default to be strict.

Sounds good, except that I don't like "casual" very much. It conveys an impression of sloppiness. Since we talk about stricness, you could have "loose", "strict", "veryStrict" (I was thinking of "paranoid", but that's probably too much ;) )

For the normal strict tearing I will set the threshold to 1e-6*p.nominal for parameters and think about some kind of appropriate error message/warning in case something bad happens after overriding.

Sounds good.

Maybe it would be enough to provide a warning that something bad could happen when a parameter is overridden with a small value or zero and tearingStrictness has not been veryStrict?

Absolutely!

[ptaeuber]

Sounds good, except that I don't like "casual" very much. It conveys an impression of sloppiness.

Hm, I thought of "casual" because this is the name I used in the context of Dynamic Tearing (strict set, casual set), where the casual set is torn with exactly that "loose strictness". So, it would probably be good to hold on to it to not introduce new phrases meaning the same as existing ones. And by the way, to allow to divide by everything, including variables changing its values, could be considered as kind of sloppy, IMHO :)

[casella]

Replying to ptaeuber:

Sounds good, except that I don't like "casual" very much. It conveys an impression of sloppiness.

Hm, I thought of "casual" because this is the name I used in the context of Dynamic Tearing (strict set, casual set), where the casual set is torn with exactly that "loose strictness". So, it would probably be good to hold on to it to not introduce new phrases meaning the same as existing ones.

Agreed.

And by the way, to allow to divide by everything, including variables changing its values, could be considered as kind of sloppy, IMHO :)

Sure :)

[ptaeuber]

In ​d379a34 --tearingStrictness=<casual|strict|veryStrict> is introduced.

In ​e36b947 a new threshold for the division by parameters is introduced. If you want to solve an equation for variable x which is multiplied by a parameter p, until now the parameter had to be at least greater than 1e-15 to allow dividing by it. Now p has to be at least 1e-6*p.nominal.

[ptaeuber]

In ​3b5a5ea a warning for overriding parameters with small values is added.

[ptaeuber]

I think the mentioned points are handled now and this ticket can be closed.