Questionable handling of initialization of a purely algebraic system with homotopy

[casella]

We are dealing with a model that we cannot share publicly, for which the backend shows weird behaviour at initialization.

The model is purely algebraic, does not contain any der() operator nor any fixed = false parameter. Hence, I would assume that the initialization problem and the simulation problem are exactly the same. However, it uses homotopy to ease the convergence of the initial equation.

The backend first complains that there are seven consistent and redundant equations in the initialization problem and removes them from the problem. This is already questionable: it is only possible to remove initial equations from the initialization problem, if they are made redundant by index reduction, but are consistent with the others. There are no initial equations in this problem, so the backend should not be allowed to remove anything.

Immediately afterwards, it complains that there are seven missing equations from the initialization problem, and then it adds fixed = true attributes to some variables, thus adding equations to the system that have nothing to do with the original problem, inevitably leading to initialization failure.

I can share the model with you privately, so you can have a look, but before doing that, I'd like you to consider some clues about what I believe are the root causes of the problem.

What I suspect is that the simplified equations introduced by homotopy at lambda = 0 make the initialization problem singular, and then the backend is trying to somehow "fix the situation". This is not correct. Singularities not involving initial equations should be reported as such, so that the modeller can fix the model; a singular model is never a good thing, except in the case of consistent initial equations made redundant by index reduction.

So, I tried to poke the backend on this point with some MWEs, from which I got a response that IMHO is highly questionable. Here they are:

package TestInitialSingular
  model M1
    Real x, y;
  equation
    sin(x) + homotopy(0.1*cos(x), 0.1) = 0;
    homotopy(sin(y), 0) = 0;
  annotation(__OpenModelica_commandLineOptions = "-d=bltdump");
  end M1;

  model M2
    Real x, y, z;
  equation
    x + y + z = 0;
    x + y + z = 0;
    x - 2*y - z = 1;
  annotation(__OpenModelica_commandLineOptions = "-d=bltdump");
end M2;
  
  model M3
    Real x, y, z;
  equation
    x + y + z = 0;
    x + homotopy(2*y,y) + z = 0;
    x - 2*y - z = 1;
  annotation(__OpenModelica_commandLineOptions = "-d=bltdump");
  end M3;

end TestInitialSingular;

Model M1 produces the following warning:

[1] 09:45:30 Translation Warning
Assuming fixed start value for the following 1 variables:
         y:VARIABLE(fixed = true ) TestInitialSingular.M1 type: Real

and then runs successfully. In my view, this is not the correct behaviour. At lambda = 0, the system's equations (no initial equation involved!) become

    sin(x) + 0.1 = 0;
    0 = 0;

This system is obviously structurally singular, because y does not appear anywhere. The root cause of this issue is an overly simplified expression in homotopy. The proper solution to this issue is reporting the singularity to the modeller, so he can try to fix the simplified expression and avoid a singular simplified system. Trying to "fix" it by adding an extra initial equation (y = y.start) is arbitrary and IMHO a very bad idea.

This is fundamentally different from adding an initial equation to this system:

model M
   Real x;
equation
  der(x) = -x + 1
end M;

This model has one state variable x, so in case the initialization problem is singular (because of lacking initial equations) it is ok to add fixed = true attributes to suitably selected variables. My point is, the backend should add fixed = true equations to fix singular initialization problems systems only if they involve state variables or fixed = false parameters, and the number of such fixed = true additional initial equations should not exceed the number of states and fixed = false parameters.

Model M2 is handled as expected by the backend, which reports:

[1] 14:51:29 Scripting Notification
unmatched equations: 1 
...
[2] 14:51:29 Translation Error
Internal error IndexReduction.pantelidesIndexReduction failed! Found empty set of continuous equations. Use -d=bltdump to get more information.

[3] 14:51:29 Symbolic Error
Model is structurally singular, error found sorting equations
  1: 0.0 = 0.0
for variables
  2: y:VARIABLE() TestInitialSingular.M2 type: Real

[4] 14:51:29 Translation Error
Internal error Transformation Module PFPlusExt index Reduction Method Pantelides failed!

I would then expect that also model M3 gives a similar outcome when analyzing the initial-lambda0 problem. What I get instead is the following warning:

[1] 14:52:14 Symbolic Warning
The linear system: 
1 : x + (-2.0) * y - z = 1.0
2 : x + y + z = 0.0
3 : x + y + z = 0.0
[
  1.0 , -1.0 , -2.0 ;
  1.0 , 1.0 , 1.0 ;
  1.0 , 1.0 , 1.0
]
  *
[
  x ;
  z ;
  y
]
  =
[
  1.0 ;
  0.0 ;
  0.0
]
 might be structurally or numerically singular for variable y since U(3,3) = 0.0. It might be hard to solve. Compilation continues anyway.

and then proceeds to initialize and simulate successfully. IMHO this is rather questionable.

Modellers introduce simplified equations in homotopy operators to ease the convergence of the initialization problem, e.g. by replacing nonlinear equations with linear ones or removing some dependencies. If the outcome of this simplification at lambda = 0 is a singular system, that's not good: they should receive feedback. Trying to solve this singular system anyway is a bad idea, because there is no control on the outcome. Initialization may succeed, but there is no guarantee that the found solution makes any sense. In this case the simulation somehow runs anyway, because the actual initial problem is linear and non-singular, so it is solved correctly no matter what the homotopy process does. However, you can very well understand that if you get a singular lambda = 0 problem and you manage to solve it somehow, it is very likely that the solution won't be the desired one, and this will later on lead to a failure on the homotopy process if seriously nonlinear equations are involved.

Singularities at lambda = 0 should be reported properly and abort the compilation process, not swept under the carpet :sweat_smile

Desired behaviour

1. The backend should never add to a (sub)system of initialization equation a number of fixed = true attributes that exceeds the number of states and parameters with fixed = false involved in that (sub)system.
2. The backend should never remove redundant equations from an initialization (sub)system, only redundant initial equations
3. Structurally singular homotopy-based initialization (sub)problems at lambda = 0 should be reported and abort the compilation process

[casella]

I discussed this with @kabdelhak. Currently there are issues in the way the lambda0 equations are handled in the old backend, stemming from the fact that the runtime has only one array for the unknown of the initialization problem, so you have to keep the same variable numbering for lambda0-initial and initial equations (I hope I understood this correctly).

He and @phannebohm are reworking and cleaning up the whole thing in the new backend. Once it's proven to work, it should be back-ported in the old one.

[kabdelhak]

I fixed it for the first model M1 in PR #10937. It now correctly reports

"[/home/kab/workspace/2_OM_Build/OMCompiler/Compiler/BackEnd/Initialization.mo:1211:5-1211:82:writable] Error: Internal error Initialization.preBalanceInitialSystem2 failed because variable y:VARIABLE() type: Real does not appear in any equation in the lambda 0 initial system and is not fixable.

As for the last model M3 we would need to use some kind of ASSC or CSE on the initial system to find these structural singularities. It is not really high on my priority list right now, working on #10920. Unfortunately the changes do not affect that model.

[phannebohm]

There is in fact the flag --initASSC which on M3 reports

[ASSC] The equation: x + y + z = 0.0
[ASSC] Gets replaced by equation: 0.0 = 0
--- Some equations have been changed, for more information please use -d=dumpASSC.---

but then keeps on compiling and simulating successfully and concludes with

Error: Initialization problem is structural singular. Please, check the initial conditions.
Warning: The initial conditions are not fully specified. For more information set -d=initialization. In OMEdit Tools->Options->Simulation->Show additional information from the initialization process, in OMNotebook call setCommandLineOptions(\"-d=initialization\").
Warning: The initial conditions are over specified. For more information set -d=initialization. In OMEdit Tools->Options->Simulation->Show additional information from the initialization process, in OMNotebook call setCommandLineOptions(\"-d=initialization\").

I guess that first error should let the model fail instead.

[casella]

I fixed it for the first model M1 in PR #10937. It now correctly reports

"[/home/kab/workspace/2_OM_Build/OMCompiler/Compiler/BackEnd/Initialization.mo:1211:5-1211:82:writable] Error: Internal error Initialization.preBalanceInitialSystem2 failed because variable y:VARIABLE() type: Real does not appear in any equation in the lambda 0 initial system and is not fixable.

Good, thanks!

[casella]

As for the last model M3 we would need to use some kind of ASSC or CSE on the initial system to find these structural singularities.

That is not necessary. My point is that if the initial system of equations is singular (also because of a numerically detected singularity), then the initialization should fail.

[casella]

There is in fact the flag --initASSC which on M3 reports

[ASSC] The equation: x + y + z = 0.0
[ASSC] Gets replaced by equation: 0.0 = 0

Obviously 0 = 0 should trigger a singularity failure.

I guess that first error should let the model fail instead.

I agree.

Additionally, the equation that gets replaced by an identity should be reported as redundant. Otherwise it is not clear how to fix the singularity.

[casella]

BTW, I also find the message

Error: Initialization problem is structural singular. Please, check the initial conditions.

highly questionable. The initialization problem can also be singular because of if initial() equations, or because of lambda = 0 equations (not initial equations)

[casella]

BTW, in case one gets a numerically singular (or close to singular, say cond > 1e10), it would be worth considering implementing the nullspace analysis I suggested long time ago in #6386. This would be really helpful to figure out the root cause of the singular behaviour, without any need of sophisticated symbolic analysis. Just basic linear algebra.

I can provide more details if needed.