GenerationOfFMUs examples in the MSL are broken because of unrecognized redundant initial equations

[Francesco Casella]

The four GenerationOfFMU models in the MSL 3.2.3:
​Modelica.Electrical.Analog.Examples.GenerationOfFMUs
​Modelica.Mechanics.Rotational.Examples.GenerationOfFMUs
​Modelica.Mechanics.Translational.Examples.GenerationOfFMUs
​Modelica.Thermal.HeatTransfer.Examples.GenerationOfFMUs
all share the same issue. They involve the initialization of high-index systems with redundant and consistent initial equations. All of them fail during analyzeInitialSystem with this kind of error

Error: The given system is mixed-determined.   [index > 3]
Please checkout the option "--maxMixedDeterminedIndex".
Error: No system for the symbolic initialization was generated

I prepared a reduced model that replicates the issue, please find it attached. It contains the following two variable declarations with fixed = true

Real directCapacity.heatCapacitor.T(quantity = "ThermodynamicTemperature", unit = "K", displayUnit = "degC", min = 0.0, start = 293.15, fixed = true, nominal = 300.0) "Temperature of element";
Real inverseCapacity.mass.T(quantity = "ThermodynamicTemperature", unit = "K", displayUnit = "degC", min = 0.0, start = 293.15, fixed = true, nominal = 300.0) "Temperature of element";

and the following chain of equalities

directCapacity.heatCapacitor.T = directCapacity.heatCapacitor.port.T;
directCapacity.heatCapacitor.port.T = directCapacity.heatFlowToTemperature.heatPort.T;
directCapacity.heatFlowToTemperature.y = directCapacity.heatFlowToTemperature.heatPort.T
directCapacity.heatFlowToTemperature.y = directCapacity.heatFlowToTemperature.p;
directCapacity.heatFlowToTemperature.p = directCapacity.T;
directCapacity.T = inverseCapacity.T;
inverseCapacity.temperatureToHeatFlow.p = inverseCapacity.T;
inverseCapacity.temperatureToHeatFlow.u = Modelica.Blocks.Interfaces.Adaptors.Functions.state1({inverseCapacity.temperatureToHeatFlow.p, inverseCapacity.temperatureToHeatFlow.u1}, time);
inverseCapacity.temperatureToHeatFlow.u = inverseCapacity.temperatureToHeatFlow.heatPort.T
inverseCapacity.temperatureToHeatFlow.heatPort.T = inverseCapacity.mass.port.T;
inverseCapacity.mass.T = inverseCapacity.mass.port.T;

which is ultimately equivalent to

directCapacity.heatCapacitor.T = inverseCapacity.mass.T

that involves two differentiated variables, thus leading to an index-2 problem with overdetermined and consistent initial equations.

We do have algorithms in OMC to deal with this case (see ​paper). The chain of equality equations is trivial and handled by alias elimination, except for the equation

inverseCapacity.temperatureToHeatFlow.u = Modelica.Blocks.Interfaces.Adaptors.Functions.state1({inverseCapacity.temperatureToHeatFlow.p, inverseCapacity.temperatureToHeatFlow.u1}, time);

The function state1() has lateInline = true, so it should be eventually eliminated once the index has been analyzed. By turning on -d=optdaedump this seems to be the case: the pre-optimization phase ends with the state1() function call still in place

7/7 (1): inverseCapacity.mass.T = Modelica.Blocks.Interfaces.Adaptors.Functions.state1({directCapacity.heatFlowToTemperature.y, directCapacity.derT}, time)   [dynamic |0|0|0|0|] 

but then, the next time the dynamic equations are printed out, it shows up (correctly) as

7/7 (1): inverseCapacity.mass.T = directCapacity.heatFlowToTemperature.y   [dynamic |0|0|0|0|] 

However, for some reason the redundant initial equation is ultimately not eliminated, leading to a mixed-determined initialization problem.

I tentatively assign this issue to @Karim.Abdelhak, since he's working on index reduction. You may co-ordinate with @lochel, who wrote the algorithm to remove the redundant initial equations.

[anonymous]

We are experiencing the same problem. Is there a way to temporary fix this problem?
Thanks.

[Francesco Casella]

Replying to anonymous:

We are experiencing the same problem. Is there a way to temporary fix this problem?
Thanks.

Of course. Just change the fixed attribute of one the two temperatures to false

[Francesco Casella]

The same issue also affects ​Modelica.Electrical.Analog.Examples.ResonanceCircuits.

Basically, all models in the MSL that try to demonstrate how to deal with index reduction across FMUs are broken.

@lochel, could you have a look at this? I gues this is very relevant w.r.t. what you are currently doing.

[anonymous]

Of course. Just change the fixed attribute of one the two temperatures to false

The problem is in the mechanical block. How can we solve this?

[Francesco Casella]

Replying to anonymous:

Of course. Just change the fixed attribute of one the two temperatures to false

The problem is in the mechanical block. How can we solve this?

Please, be more specific. What is "the mechanical block"? Do you mean you have the same problem as in Modelica.Mechanics.Rotational.Examples.GenerationOfFMUs?

Ideally, it would be nice if you could attach a test case to this ticket that demonstrates your problem.

[Francesco Casella]

OMC 2.0.0 should run all the MSL correctly, so I'm marking this as a blocker for that version.

[Francesco Casella]

Try to get 100% MSL 3.2.3 to work correctly for 1.14.0

[Karim Adbdelhak]

As lennart commented on ​PR3085 my changes are incorrect because the introduced artificial equations should only be connected to inInitVarIndices (states, discrete states, and parameters having fixed=false)

Here inInitVarIndices always resolves to be empty, i guess it is not correctly determined. looking into it!

[Lennart Ochel]

This is currently under investigation.

[Francesco Casella]

@Karim.Abdelhak, @lochel, can you please update on the status of this issue? Please note that PR3085 is on a dead track, because the OMCompiler repository was archived. So you need to resurrect it on the OpenModelica repository if you want to continue with it.

[Lennart Ochel]

After today’s discussion with Karim, we think it might be related to the pre-optimization module resolveLoops.

[Karim Adbdelhak]

With the old frontend everything works when using --preOptModules-=evaluateParameters,simplifyIfEquations, we did not find any way to make it work with -d=newInst. As it seems there is a difference in detecting structural parameters in bindings of other parameters and handling of protected parameters.

Related discussion for that can be found at ​PR318.

[Karim Adbdelhak]

The old frontend behaves differently from the new frontend here, thats why i added @perost.

There is some difference in the evaluation of parameters.

[Francesco Casella]

@Karim I'm trying to make sense of this ticket, but I need some help from you.

The conclusion of my initial analysis was that in all these cases the problem was that the algorithm that detects and removes redundant initial equations caused by index reduction was not working properly, as demonstrated by the attached test case. Can you confirm this is still correct?

That said, I see that the OF returns:

  parameter Boolean directCapacity.heatFlowToTemperature.use_pder = true 
    "Use output for 1st derivative of potential" annotation(Evaluate = true);
equation
  directCapacity.heatFlowToTemperature.y1 =
   if directCapacity.heatFlowToTemperature.use_pder then   
     der(directCapacity.heatFlowToTemperature.y) 
   else 
     0.0;

while the NF returns

  final parameter Boolean directCapacity.heatFlowToTemperature.use_pder = true "Use output for 1st derivative of potential" annotation(Evaluate = true);
equation
  directCapacity.heatFlowToTemperature.y1 =
   der(directCapacity.heatFlowToTemperature.y);

It seems to me that the NF is doing the right thing: use_pder has annotation(Evaluate = true), so the NF exploits that to simplify away the conditional expression.

By the way, use_pder is actually a structural parameter, even though the frontend may not recognize this, because if it is set to false, there are no more instances of der(directCapacity.heatFlowToTemperature.y) in the model, so the order of the system would change. That's the reason why annotation(Evaluate = true) was put there in the first place, to make sure the conditional expression doesn't make it to the code generation phase even if the compiler does not recognize the structural attribute of use_pder.

Why keeping the conditional expression actually helps the backend to figure out the redundant initial equation, I have no idea, but for sure the NF is doing the right thing. If the current backend is more comfortable with the conditional expression (which is potentially destructive from a structural point of view), that's a problem of the backend.

BTW, I also tried to run the attached test case with -d=newInst --preOptModules-=evaluateParameters,simplifyIfEquations and the result was that after two minutes OMEdit was still in the process of elaborating something, at which point I just killed it. That puzzles me even more.

[Lennart Ochel]

After spending a lot of time on this issue, we came up with a minimal model showing the problem:

model M
  function state1
    input Real u[2];
    input Real dummy;
    output Real s;
  algorithm
    s := u[1];
  annotation (derivative(noDerivative=u) = state1der1, InlineAfterIndexReduction=true);
  end state1;

  function state1der1
    input Real u[2];
    input Real dummy;
    input Real dummy_der;
    output Real sder1;
  algorithm
    sder1 := u[2];
  annotation (InlineAfterIndexReduction=true);
  end state1der1;

  parameter Real C1 = 10;
  parameter Real C2 = 20;
  Real c1v, c1i, c2v, c2i, i;
equation
  c2v = state1({c1v, der(c1v)},time);
  i = sin(time);
  i = c1i + c2i;
  der(c1v) = 1/C1 * c1i;
  der(c2v) = 1/C2 * c2i;
end M;

The problem is related to the functions and the InlineAfterIndexReduction annotation.

[Karim Adbdelhak]

In addition to what @lochel said i will point out the following:

If we match the system above without inlining, equation 1

c2v = state1({c1v, der(c1v)},time);

gets marked to be (implicitly) solvable by der(c1v). This is not true and would be clear if we inlined the function beforehand, because it just resolves to be

c2v = c1v;

InlineAfterIndexReduction=true prevents that though. Because of this misinformation the compiler does not find the MSSS in this case, which would just be equation 1.

For index reduction, on the other hand, we need to have it not inlined, such that the compiler knows that the derivative of state1 will be state1der1. In this case it would work inlined, but in the proposed original Models it would not.

It seems like we need two different versions of this equation for matching and deriving. A better name for the attribute InlineAfterIndexReduction would be something like DoNotInlineForDifferentiation or something similar i guess. We have to come up with a consistent solution for this. One solution would be to prevent implicit solving of functions in general, but i don't think that would be the right way to go.

[Francesco Casella]

Replying to Karim.Abdelhak:

One solution would be to prevent implicit solving of functions in general, but i don't think that would be the right way to go.

What about preventing the implicit solutions of functions that have the InlineAfterIndexReduction annotation?

I scanned the whole MSL 3.2.3, there are just 9 such functions:

• 5 in Modelica.Blocks.Interfaces, which are the index-2 and index-3 versions of the function used here
• 3 in Modelica.Mechanics.MultiBody, which are resolve1, resolve2, and resolveRelative
• 1 in Modelica.Mechanics.Rotational, function Move_w (not sure why there is no corresponding function in Modelia.Mechanics.Translational, probably an oversight).

As I understand, the use pattern is similar, so I guess it could work.

[Karim Adbdelhak]

Replying to casella:

What about preventing the implicit solutions of functions that have the InlineAfterIndexReduction annotation?

Well i guess it is a common case for functions which have been marked with InlineAfterIndexReduction to contain more dependencies than they actually have to correctly construct the derivative. This still feels like a hack since there may be models where you would want to have it be implicitely solved for a variable. Maybe an inlined version should be created for matching, without putting it in the actual system.

I am not quite sure, this seems like a corner case i guess, but other libraries may use this quite often. If i have time next monday i will try it.

[Karim Adbdelhak]

Inlining function calls with InlineAfterIndexReduction before and only for matching does the trick. Inlining these functions is now done twice (before and after index recuction). Since it does not occure that often it should not make much of a difference regarding compilation speed, saving and passing around the result could even be slower.

Changes can be seen in ​PR326.

[Karim Adbdelhak]

All models can be compiled and simulated now, although the electrical and thermal model seem to have an overdetermination which apparently can't be checked by the backend.

@casella Could you check the results and verify them?

[Karim Adbdelhak]

Looking at the latest test result changes from ​2019-07-15 we can see that most of it is resolved, but some are not verified jet. Mixed results regarding simulation time, but mostly positive.

It also seems like some models containing Joints broke, maybe we need to investigate that.

[Francesco Casella]

Replying to Karim.Abdelhak:

All models can be compiled and simulated now, although the electrical and thermal model seem to have an overdetermination which apparently can't be checked by the backend.

Where can I see that?

[Karim Adbdelhak]

Replying to casella:

Replying to Karim.Abdelhak:

All models can be compiled and simulated now, although the electrical and thermal model seem to have an overdetermination which apparently can't be checked by the backend.

Where can I see that?

I ran follwing .mos file:

loadModel(Modelica,{"trunk"});getErrorString();
simulate(Modelica.Thermal.HeatTransfer.Examples.GenerationOfFMUs);getErrorString();

and got the warning

"Warning: It was not possible to determine if the initialization problem is consistent, because of not evaluable parameters/start values during compile time: directCapacity.heatFlowToTemperature.y = $START.directCapacity.heatFlowToTemperature.y ($START.inverseCapacity.mass.T = $START.directCapacity.heatFlowToTemperature.y)
Warning: The initial conditions are over specified. For more information set -d=initialization. In OMEdit Tools->Options->Simulation->OMCFlags, in OMNotebook call setCommandLineOptions("-d=initialization").
"

[Francesco Casella]

Replying to Karim.Abdelhak:

I ran following .mos file:

loadModel(Modelica,{"trunk"});getErrorString();
simulate(Modelica.Thermal.HeatTransfer.Examples.GenerationOfFMUs);getErrorString();

and got the warning

"Warning: It was not possible to determine if the initialization problem is consistent, because of not evaluable parameters/start values during compile time: directCapacity.heatFlowToTemperature.y = $START.directCapacity.heatFlowToTemperature.y ($START.inverseCapacity.mass.T = $START.directCapacity.heatFlowToTemperature.y)
Warning: The initial conditions are over specified. For more information set -d=initialization. In OMEdit Tools->Options->Simulation->OMCFlags, in OMNotebook call setCommandLineOptions("-d=initialization").
"

Well, as explained in the original description of this ticket, the second warning is indeed correct. It seems that the backend algorithm dealing with redundant initial equations can now handle that, given that the simulation runs and passes the verification stage. In fact if you run the model with -d=initialization you get the following message:

[2] 13:57:19 Translation Warning
The initial conditions are over specified. The following 1 initial equations are redundant, so they are removed from the initialization sytem:
         directCapacity.heatFlowToTemperature.y = $START.directCapacity.heatFlowToTemperature.y.

which is precisely the expected behaviour. I'm not sure right now about the first warning, but it doesn't worry me too much at this moment.

[Francesco Casella]

Replying to Karim.Abdelhak:

Looking at the latest test result changes from ​2019-07-15 we can see that most of it is resolved, but some are not verified jet.

Most models not pass verification don't do so because the reference results are missing. The only one that actually fails is ​Modelica.Mechanics.Rotational.Examples.GenerationOfFMUs (in all versions of MSL). I checked, and it seems that the OMC-generated solutions barely touches the tolerance tube at some points in the simulation. I would classify this as a false negative, we'll see later on how to deal with these cases, there are quite many in the MSL. Don't worry about that.

Mixed results regarding simulation time, but mostly positive.

This measurement is not reliable at all, due to some reasons (probably because of memory access issues when memory simulations run in parallel), simulation times often get doubled or halved without any obvious reasons (e.g., when OMEdit-only changes were pushed). Don't worry about that either.

It also seems like some models containing Joints broke, maybe we need to investigate that.

Yeah, that's the only thing worth checking. In fact, at the end of the day we only have two cases:

1. ModelicaTest.MultiBody.Joints.RevolutePlanarLoopConstraint fails with the C runtime due to some ​singularity at initialization
2. ModelicaTest.MultiBody.InitializationConversion.Joints fails with the C++ runtime due to ​some issue at initialization.

Do you want to have a look at this yourself?

[Karim Adbdelhak]

Replying to casella:

Replying to Karim.Abdelhak:

It also seems like some models containing Joints broke, maybe we need to investigate that.

Yeah, that's the only thing worth checking. In fact, at the end of the day we only have two cases:

1. ModelicaTest.MultiBody.Joints.RevolutePlanarLoopConstraint fails with the C runtime due to some ​singularity at initialization
2. ModelicaTest.MultiBody.InitializationConversion.Joints fails with the C++ runtime due to ​some issue at initialization.

Do you want to have a look at this yourself?

This week i will be occupied with grading exams and only have very little time for anything else and following two weeks i will be on vacation. If you got time i would be glad about some analytics from your point of view so i can pick up on that when i am back.

[Francesco Casella]

OK. I'm going to close this ticket because the original issue is now solved for good.

For some analysis of the remaining issues, see #5590.