More issues with dynamic state selection and Jacobians

[casella]

This issue is related to #5459 but I'm opening a separate ticket here because that ticket is too cluttered now.

Consider the attached test case, which is part of the FlexiCaL plant model that I'd like to present at the next Annual Meeting and to use as showcase for OMC in industrial-grade power generation systems, which is possible because the model is 100% public.

Dymola produces the intended fixed set of 151 states

Statically selected continuous time states
fWPH_12_1.flowDRAIN.h[2]
fWPH_12_1.flowDRAIN.h[3]
fWPH_12_1.flowDRAIN.h[4]
fWPH_12_1.flowDRAIN.h[5]
fWPH_12_1.flowDRAIN.h[6]
fWPH_12_1.flowDRAIN.p
fWPH_12_1.flowFW_cond.htilde[]
fWPH_12_1.flowFW_cond.p
fWPH_12_1.flowFW_des.htilde[]
fWPH_12_1.flowFW_des.p
fWPH_12_1.flowFW_sub.htilde[]
fWPH_12_1.flowFW_sub.p
fWPH_12_1.flowSTEAM.h[2]
fWPH_12_1.flowSTEAM.h[3]
fWPH_12_1.flowSTEAM.h[4]
fWPH_12_1.flowSTEAM.h[5]
fWPH_12_1.flowSTEAM.h[6]
fWPH_12_1.Level_PID.I
fWPH_12_1.metalTubeFV_cond.Tvol[]
fWPH_12_1.metalTubeFV_des.Tvol[]
fWPH_12_1.metalTubeFV_sub.Tvol[]
fWPH_12_1.nusseltCondenser.hl
fWPH_12_1.nusseltCondenser.hv
fWPH_12_1.nusseltCondenser.p
fWPH_12_1.nusseltCondenser.zl
fWPH_12_1.ValveDynamic.y
fWPH_13_1.flowDRAIN.h[2]
fWPH_13_1.flowDRAIN.h[3]
fWPH_13_1.flowDRAIN.h[4]
fWPH_13_1.flowDRAIN.h[5]
fWPH_13_1.flowDRAIN.h[6]
fWPH_13_1.flowDRAIN.p
fWPH_13_1.flowFW_cond.htilde[]
fWPH_13_1.flowFW_cond.p
fWPH_13_1.flowFW_des.htilde[]
fWPH_13_1.flowFW_des.p
fWPH_13_1.flowFW_sub.htilde[]
fWPH_13_1.flowFW_sub.p
fWPH_13_1.flowSTEAM.h[2]
fWPH_13_1.flowSTEAM.h[3]
fWPH_13_1.flowSTEAM.h[4]
fWPH_13_1.flowSTEAM.h[5]
fWPH_13_1.flowSTEAM.h[6]
fWPH_13_1.Level_PID.I
fWPH_13_1.metalTubeFV_cond.Tvol[]
fWPH_13_1.metalTubeFV_des.Tvol[]
fWPH_13_1.metalTubeFV_sub.Tvol[]
fWPH_13_1.nusseltCondenser.hl
fWPH_13_1.nusseltCondenser.hv
fWPH_13_1.nusseltCondenser.p
fWPH_13_1.nusseltCondenser.zl
fWPH_13_1.ValveDynamic.y
fWPH_14_1.flowDRAIN.h[2]
fWPH_14_1.flowDRAIN.h[3]
fWPH_14_1.flowDRAIN.h[4]
fWPH_14_1.flowDRAIN.h[5]
fWPH_14_1.flowDRAIN.h[6]
fWPH_14_1.flowDRAIN.p
fWPH_14_1.flowFW_cond.htilde[]
fWPH_14_1.flowFW_cond.p
fWPH_14_1.flowFW_des.htilde[]
fWPH_14_1.flowFW_des.p
fWPH_14_1.flowFW_sub.htilde[]
fWPH_14_1.flowFW_sub.p
fWPH_14_1.flowSTEAM.h[2]
fWPH_14_1.flowSTEAM.h[3]
fWPH_14_1.flowSTEAM.h[4]
fWPH_14_1.flowSTEAM.h[5]
fWPH_14_1.flowSTEAM.h[6]
fWPH_14_1.Level_PID.I
fWPH_14_1.metalTubeFV_cond.Tvol[]
fWPH_14_1.metalTubeFV_des.Tvol[]
fWPH_14_1.metalTubeFV_sub.Tvol[]
fWPH_14_1.nusseltCondenser.hl
fWPH_14_1.nusseltCondenser.hv
fWPH_14_1.nusseltCondenser.p
fWPH_14_1.nusseltCondenser.zl
fWPH_14_1.ValveDynamic.y
IP_turbine3.phi

and runs the simulation in about 2 s.

OMC instead produces the following dynamic state selection:

 * Number of states: 151 ($STATESET1.x[6],$STATESET1.x[5],$STATESET1.x[4],
 $STATESET1.x[3],$STATESET1.x[2],$STATESET1.x[1],
constantSpeed.phi,
fWPH_12_1.ValveDynamic.y,
fWPH_12_1.Level_PID.I,
fWPH_12_1.flowFW_des.p,
fWPH_12_1.flowFW_des.htilde[1],
fWPH_12_1.flowFW_des.htilde[2],
fWPH_12_1.flowFW_des.htilde[3],
fWPH_12_1.flowFW_des.htilde[4],
fWPH_12_1.flowFW_des.htilde[5],
fWPH_12_1.flowFW_cond.p,
fWPH_12_1.flowFW_cond.htilde[1],
fWPH_12_1.flowFW_cond.htilde[2],
fWPH_12_1.flowFW_cond.htilde[3],
fWPH_12_1.flowFW_cond.htilde[4],
fWPH_12_1.flowFW_cond.htilde[5],
fWPH_12_1.flowFW_sub.p,
fWPH_12_1.flowFW_sub.htilde[1],
fWPH_12_1.flowFW_sub.htilde[2],
fWPH_12_1.flowFW_sub.htilde[3],
fWPH_12_1.flowFW_sub.htilde[4],
fWPH_12_1.flowFW_sub.htilde[5],
fWPH_12_1.flowDRAIN.p,
fWPH_12_1.flowDRAIN.htilde[1],
fWPH_12_1.flowDRAIN.htilde[2],
fWPH_12_1.flowDRAIN.htilde[3],
fWPH_12_1.flowDRAIN.htilde[4],
fWPH_12_1.flowDRAIN.htilde[5],
fWPH_12_1.flowSTEAM.htilde[1],
fWPH_12_1.flowSTEAM.htilde[2],
fWPH_12_1.flowSTEAM.htilde[3],
fWPH_12_1.flowSTEAM.htilde[4],
fWPH_12_1.flowSTEAM.htilde[5],
fWPH_12_1.metalTubeFV_des.Tvol[1],
fWPH_12_1.metalTubeFV_des.Tvol[2],
fWPH_12_1.metalTubeFV_des.Tvol[3],
fWPH_12_1.metalTubeFV_des.Tvol[4],
fWPH_12_1.metalTubeFV_des.Tvol[5],
fWPH_12_1.metalTubeFV_cond.Tvol[1],
fWPH_12_1.metalTubeFV_cond.Tvol[2],
fWPH_12_1.metalTubeFV_cond.Tvol[3],
fWPH_12_1.metalTubeFV_cond.Tvol[4],
fWPH_12_1.metalTubeFV_cond.Tvol[5],
fWPH_12_1.metalTubeFV_sub.Tvol[1],
fWPH_12_1.metalTubeFV_sub.Tvol[2],
fWPH_12_1.metalTubeFV_sub.Tvol[3],
fWPH_12_1.metalTubeFV_sub.Tvol[4],
fWPH_12_1.metalTubeFV_sub.Tvol[5],
fWPH_12_1.nusseltCondenser.zl,
fWPH_12_1.nusseltCondenser.p,
fWPH_12_1.nusseltCondenser.hv,
fWPH_12_1.nusseltCondenser.hl,
fWPH_13_1.ValveDynamic.y,
fWPH_13_1.Level_PID.I,
fWPH_13_1.flowFW_des.p,
fWPH_13_1.flowFW_des.htilde[1],
fWPH_13_1.flowFW_des.htilde[2],
fWPH_13_1.flowFW_des.htilde[3],
fWPH_13_1.flowFW_des.htilde[4],
fWPH_13_1.flowFW_des.htilde[5],
fWPH_13_1.flowFW_cond.p,
fWPH_13_1.flowFW_cond.htilde[1],
fWPH_13_1.flowFW_cond.htilde[2],
fWPH_13_1.flowFW_cond.htilde[3],
fWPH_13_1.flowFW_cond.htilde[4],
fWPH_13_1.flowFW_cond.htilde[5],
fWPH_13_1.flowFW_sub.p,
fWPH_13_1.flowFW_sub.htilde[1],
fWPH_13_1.flowFW_sub.htilde[2],
fWPH_13_1.flowFW_sub.htilde[3],
fWPH_13_1.flowFW_sub.htilde[4],
fWPH_13_1.flowFW_sub.htilde[5],
fWPH_13_1.flowDRAIN.p,
fWPH_13_1.flowDRAIN.htilde[1],
fWPH_13_1.flowDRAIN.htilde[2],
fWPH_13_1.flowDRAIN.htilde[3],
fWPH_13_1.flowDRAIN.htilde[4],
fWPH_13_1.flowDRAIN.htilde[5],
fWPH_13_1.metalTubeFV_des.Tvol[1],
fWPH_13_1.metalTubeFV_des.Tvol[2],
fWPH_13_1.metalTubeFV_des.Tvol[3],
fWPH_13_1.metalTubeFV_des.Tvol[4],
fWPH_13_1.metalTubeFV_des.Tvol[5],
fWPH_13_1.metalTubeFV_cond.Tvol[1],
fWPH_13_1.metalTubeFV_cond.Tvol[2],
fWPH_13_1.metalTubeFV_cond.Tvol[3],
fWPH_13_1.metalTubeFV_cond.Tvol[4],
fWPH_13_1.metalTubeFV_cond.Tvol[5],
fWPH_13_1.metalTubeFV_sub.Tvol[1],
fWPH_13_1.metalTubeFV_sub.Tvol[2],
fWPH_13_1.metalTubeFV_sub.Tvol[3],
fWPH_13_1.metalTubeFV_sub.Tvol[4],
fWPH_13_1.metalTubeFV_sub.Tvol[5],
fWPH_13_1.nusseltCondenser.zl,
fWPH_13_1.nusseltCondenser.hv,
fWPH_13_1.nusseltCondenser.hl,
fWPH_14_1.ValveDynamic.y,
fWPH_14_1.Level_PID.I,
fWPH_14_1.flowFW_des.p,
fWPH_14_1.flowFW_des.htilde[1],
fWPH_14_1.flowFW_des.htilde[2],
fWPH_14_1.flowFW_des.htilde[3],
fWPH_14_1.flowFW_des.htilde[4],
fWPH_14_1.flowFW_des.htilde[5],
fWPH_14_1.flowFW_cond.p,
fWPH_14_1.flowFW_cond.htilde[1],
fWPH_14_1.flowFW_cond.htilde[2],
fWPH_14_1.flowFW_cond.htilde[3],
fWPH_14_1.flowFW_cond.htilde[4],
fWPH_14_1.flowFW_cond.htilde[5],
fWPH_14_1.flowFW_sub.p,
fWPH_14_1.flowFW_sub.htilde[1],
fWPH_14_1.flowFW_sub.htilde[2],
fWPH_14_1.flowFW_sub.htilde[3],
fWPH_14_1.flowFW_sub.htilde[4],
fWPH_14_1.flowFW_sub.htilde[5],
fWPH_14_1.flowDRAIN.p,
fWPH_14_1.flowDRAIN.htilde[1],
fWPH_14_1.flowDRAIN.htilde[2],
fWPH_14_1.flowDRAIN.htilde[3],
fWPH_14_1.flowDRAIN.htilde[4],
fWPH_14_1.flowDRAIN.htilde[5],
fWPH_14_1.flowSTEAM.htilde[1],
fWPH_14_1.flowSTEAM.htilde[2],
fWPH_14_1.flowSTEAM.htilde[3],
fWPH_14_1.flowSTEAM.htilde[4],
fWPH_14_1.flowSTEAM.htilde[5],
fWPH_14_1.metalTubeFV_des.Tvol[1],
fWPH_14_1.metalTubeFV_des.Tvol[2],
fWPH_14_1.metalTubeFV_des.Tvol[3],
fWPH_14_1.metalTubeFV_des.Tvol[4],
fWPH_14_1.metalTubeFV_des.Tvol[5],
fWPH_14_1.metalTubeFV_cond.Tvol[1],
fWPH_14_1.metalTubeFV_cond.Tvol[2],
fWPH_14_1.metalTubeFV_cond.Tvol[3],
fWPH_14_1.metalTubeFV_cond.Tvol[4],
fWPH_14_1.metalTubeFV_cond.Tvol[5],
fWPH_14_1.metalTubeFV_sub.Tvol[1],
fWPH_14_1.metalTubeFV_sub.Tvol[2],
fWPH_14_1.metalTubeFV_sub.Tvol[3],
fWPH_14_1.metalTubeFV_sub.Tvol[4],
fWPH_14_1.metalTubeFV_sub.Tvol[5],
fWPH_14_1.nusseltCondenser.zl,
fWPH_14_1.nusseltCondenser.p,
fWPH_14_1.nusseltCondenser.hv,
fWPH_14_1.nusseltCondenser.hl)

with 6 dynamically chosen states.

This has two very bad consequences: a very long compilation time and a very long simulation time, around 30 seconds. I suspect both are due to the unnecessary need of solving some large algebraic loops due to the bad choice of states, since the nls.c file is 12 MB, at least 6X bigger than all the other ones, and it is by far the longest to compile.

The lenghty simulation may also be traced back to the very long time taken to compute the numerical Jacobians, around 8 seconds. All results are obtained with an i7-8550 Intel CPU.

During the simulation, we get the error message

Warning: maximal number of iteration reached but no root found

about 120 times, which also suggests something wrong with the state selection is unnecessarily

I'm not sure if coloured Jacobians are not available in case of dynamic state selection, which could partly explain this very long time, since the system has 151 states. I ask for confirmation of this possible explanation.

In any case, it is well-known from other test cases (such as those reported in #5459) that having a dynamic state selection instead of the right static one often implies a huge performance degradation.

@Karim, we should be able to get the same static selection that Dymola gets (or something equivalent in terms of aliases), which is also the one I expect from a modelling point of view.

Getting this model to work properly has strategic value for me and for OMC in general, so I'd be glad if you could have a look at it. I think the issues involved are still the same that plagued #5459.

Thanks!

[casella]

Test script

[casella]

Save Total model

[Karim.Abdelhak]

From first glance i can see that we select variable 1-5 from e.g. fWPH_12_1.flowDRAIN.htilde whereas dymola selects 2-6. This pattern seems to repeat over multiple of those arrays, this could be the source of the problem. Do you have any intel on why these selections are mandatory?

I will try to have a deeper look at it this week!

[casella]

Replying to Karim.Abdelhak:

From first glance i can see that we select variable 1-5 from e.g. fWPH_12_1.flowDRAIN.htilde whereas dymola selects 2-6.

I don't think that is the point. The models are structured in this way

model M
  parameter Integer N;
  SI.SpecificEnthalpy h[N];
  SI.SpecificEnthalpy htilde[N-1];
  ...
equation
 ...
 htilde = h[2:N];
 ...
 for i in 1:N-1 loop
   // equations involving der(htilde)
 end for;
 ...
end M;

So, htilde[1:5] and h[2:6] are just aliases.

[Karim.Abdelhak]

I think i found the root of the problem. We have following subset of constraint equations during index reduction:

IP_turbine1.w = homotopy(IP_turbine1.Kt * ParArc1 * sqrt(sourcePressure.p0 * IP_turbine1.rho) * (1.0 - (1.0 / IP_turbine1.PR) ^ 2.0) / ((1.0 - (1.0 / IP_turbine1.PR) ^ 2.0) ^ 2.0 + 0.0001) ^ 0.25, ParArc1 * IP_turbine1.wnom / IP_turbine1.pnom * sourcePressure.p0)
IP_turbine2.w = homotopy(IP_turbine2.Kt * ParArc2 * sqrt(IP_turbine2.pin * IP_turbine2.rho) * (1.0 - (1.0 / IP_turbine2.PR) ^ 2.0) / ((1.0 - (1.0 / IP_turbine2.PR) ^ 2.0) ^ 2.0 + 0.0001) ^ 0.25, ParArc2 * IP_turbine2.wnom / IP_turbine2.pnom * IP_turbine2.pin)
w_Bleeding1 = homotopy(valveVap1.theta_act * valveVap1.Av * valveVap1.Y * sqrt(valveVap1.rho) * v     alveVap1.p * valveVap1.xs / ((valveVap1.p * valveVap1.xs) ^ 2.0 + (valveVap1.b * valveVap1.dpnom) ^ 2.0) ^ 0.25, valveVap1.theta_act / valveVap1.thetanom * valveVap1.wnom / valveVap1.dpnom * (IP_turbine2.pin - valveVap1.outlet.p))
IP_turbine2.w + w_Bleeding1 - IP_turbine1.w = 0.0

Since for the homotopy we only consider the actual expression and not the simplified, this can be reduced to:

IP_turbine1.w = IP_turbine1.Kt * ParArc1 * sqrt(sourcePressure.p0 * IP_turbine1.rho) * (1.0 - (1.0 / IP_turbine1.PR) ^ 2.0) / ((1.0 - (1.0 / IP_turbine1.PR) ^ 2.0) ^ 2.0 + 0.0001) ^ 0.25
IP_turbine2.w = IP_turbine2.Kt * ParArc2 * sqrt(IP_turbine2.pin * IP_turbine2.rho) * (1.0 - (1.0 / IP_turbine2.PR) ^ 2.0) / ((1.0 - (1.0 / IP_turbine2.PR) ^ 2.0) ^ 2.0 + 0.0001) ^ 0.25
w_Bleeding1 = valveVap1.theta_act * valveVap1.Av * valveVap1.Y * sqrt(valveVap1.rho) * valveVap1.p * valveVap1.xs / ((valveVap1.p * valveVap1.xs) ^ 2.0 + (valveVap1.b * valveVap1.dpnom) ^ 2.0) ^ 0.25
IP_turbine2.w + w_Bleeding1 - IP_turbine1.w = 0.0

The linear appearing state candidates for these equations are:

IP_turbine1.w
IP_turbine2.w
w_Bleeding1

As we can see there are only three linear candidates in four equations, therefore the compiler concludes that dynamic state selection is needed. Further state candidates appearing nonlinear are:

IP_turbine1.rho
IP_turbine1.PR
IP_turbine2.pin
IP_turbine2.rho
IP_turbine2.PR
valveVap1.theta_act
valveVap1.Y
valveVap1.rho
valveVap1.p
valveVap1.xs

From what i can see none of these are chosen by dymola, but maybe there are some alias eliminations i am not aware of. If that is not the case, dymola does not detect these equations to be part of the index reduction problem at all.

@Francesco Do you think alias elimination is involved here, or do you think these equations should not be part of the index reduction problem at all? From what i can see it may be the case, that we might still analyze the simplified expression at some point during causalization such that it leads to this error. Some variables like IP_turbine1.wnom only appear in the simplified expression and should not be considered for matching.

[casella]

Replying to Karim.Abdelhak:

IP_turbine1.w = IP_turbine1.Kt * ParArc1 * sqrt(sourcePressure.p0 * IP_turbine1.rho) * (1.0 - (1.0 / IP_turbine1.PR) ^ 2.0) / ((1.0 - (1.0 / IP_turbine1.PR) ^ 2.0) ^ 2.0 + 0.0001) ^ 0.25
IP_turbine2.w = IP_turbine2.Kt * ParArc2 * sqrt(IP_turbine2.pin * IP_turbine2.rho) * (1.0 - (1.0 / IP_turbine2.PR) ^ 2.0) / ((1.0 - (1.0 / IP_turbine2.PR) ^ 2.0) ^ 2.0 + 0.0001) ^ 0.25
w_Bleeding1 = valveVap1.theta_act * valveVap1.Av * valveVap1.Y * sqrt(valveVap1.rho) * valveVap1.p * valveVap1.xs / ((valveVap1.p * valveVap1.xs) ^ 2.0 + (valveVap1.b * valveVap1.dpnom) ^ 2.0) ^ 0.25
IP_turbine2.w + w_Bleeding1 - IP_turbine1.w = 0.0

I still have trouble understanding why this set of equations is considered as a set of potential state constraints in the first place, and how those variables could be considered as potential states.

As far as I understand, according to Pantelides' original formulation, this set is not an MSSS; with reference to equations (9) and (10) of his famous paper, k = 4 and l = 13, so the sufficient condition for singularity (i.e. structural singularity) does not hold.

Of course if you rule out all the variables appearing nonlinearly, just in case, than l = 3 and you have potential numerical singularity, not structural singularity.

As a modeller, I am perfectly fine with that. For example, if IP_turbine1.rho ever became zero, the whole thing would come crashing down from a mathematical point of view. However, as a modeller, I know that the density is always going to be positive, so I don't care about that.

I somehow have the feeling that the implementation in OMC is trying to be too cautious to avoid potential singularities that will never materialize, at the cost of completely useless dynamic state selections that have a huge cost in terms of numerical performance.

Do I miss something here?

[casella]

I just wonder: some of the components in this model have (appropriate) StateSelect.prefer attributes, set (if I am not mistaken) on variables that actually appear under derivative operator. Some other components have differentiated variables that should actually be picked as states, but they do not have an explicit StateSelect.prefer attribute.

A fixed choice of states including the ones with StateSelect.prefer and the ones which appear under derivative operator with no StateSelect indication would be just fine, and is in fact what Dymola does. Could it be that the partial indication of StateSelect somehow confuses the algorithm?

[casella]

@Francesco Do you think alias elimination is involved here, or do you think these equations should not be part of the index reduction problem at all?

Most definitely they shouldn't, see above.

From what i can see it may be the case, that we might still analyze the simplified expression at some point during causalization such that it leads to this error.

OMG. All instances of homotopy should be replaced by the actual expression before even starting to tackle causalization. Isn't this the case already?

Homotopy is exclusively relevant for the handling of the initialization problem, which I understand is considered after causalization, when you know who the states are. The backend should restart from scratch at this point, by joining the set of the original equations (including homotopy) with the set of initial equations. Then, a further version of the problem for the lambda = 0 case should be generated, in which all the structural simplifications enabled by that condition are exploited, to facilitate the convergence of the nonlinear solvers at the first step when lambda = 0.

Do I miss something?

Some variables like IP_turbine1.wnom only appear in the simplified expression and should not be considered for matching.

IP_turbine1.wnom is a parameter, not a variable. If I am not mistaken, in this case the simplified equations involve a strict subset of the unknowns of the actual problem.

Hope this helps :)

[Karim.Abdelhak]

Replying to casella:

Replying to Karim.Abdelhak:

IP_turbine1.w = IP_turbine1.Kt * ParArc1 * sqrt(sourcePressure.p0 * IP_turbine1.rho) * (1.0 - (1.0 / IP_turbine1.PR) ^ 2.0) / ((1.0 - (1.0 / IP_turbine1.PR) ^ 2.0) ^ 2.0 + 0.0001) ^ 0.25
IP_turbine2.w = IP_turbine2.Kt * ParArc2 * sqrt(IP_turbine2.pin * IP_turbine2.rho) * (1.0 - (1.0 / IP_turbine2.PR) ^ 2.0) / ((1.0 - (1.0 / IP_turbine2.PR) ^ 2.0) ^ 2.0 + 0.0001) ^ 0.25
w_Bleeding1 = valveVap1.theta_act * valveVap1.Av * valveVap1.Y * sqrt(valveVap1.rho) * valveVap1.p * valveVap1.xs / ((valveVap1.p * valveVap1.xs) ^ 2.0 + (valveVap1.b * valveVap1.dpnom) ^ 2.0) ^ 0.25
IP_turbine2.w + w_Bleeding1 - IP_turbine1.w = 0.0

I still have trouble understanding why this set of equations is considered as a set of potential state constraints in the first place, and how those variables could be considered as potential states.

As far as I understand, according to Pantelides' original formulation, this set is not an MSSS; with reference to equations (9) and (10) of his famous paper, k = 4 and l = 13, so the sufficient condition for singularity (i.e. structural singularity) does not hold.

Of course if you rule out all the variables appearing nonlinearly, just in case, than l = 3 and you have potential numerical singularity, not structural singularity.

The equations i posted are only a subset of the full MSSS. The full thing is too large therefore i only post the equation indices to showcase this:

1233, 1175, 1177, 1262, 1171, 56, 29, 55, 30, 36, 31, 32, 48, 2, 49, 1172, 1182, 1191, 1200, 1209, 38, 43, 45, 47, 41, 42, 40, 39, 46, 1180, 1183, 1218, 1219, 1239, 1216, 1217, 1224, 1214, 1215, 1189, 1192, 1245, 1198, 1201, 1251, 1207, 1210, 1257, 1893, 1840, 1866, 1867, 1887, 1864, 1865, 1837, 1872, 1862, 1863, 1839, 1910, 18     19, 3, 65, 51, 50, 1820, 1825, 1830, 1848, 1857, 66, 59, 60, 57, 61, 58, 62, 64, 67, 1823, 1881, 1828, 1831, 1846, 1849, 1899, 1855, 1858, 1905

From these equation we try to find a matching to linear appearing states such that they can be chosen statically as dummy states, all remaining variables will be continuous states. The subset i posted cannot be solved for enough linear appearing variables and the compiler will therefore decide to chose the remaining dummy states dynamically (and therefore also the continuous ones dynamically). The rest can be fully chosen statically.

This static selection has proven to be too strict in the past, therefore i adapted it such that it forces all StateSelect.never states to be matched here, hence to be dummy states, even if they only appear nonlinear.

As a modeller, I am perfectly fine with that. For example, if IP_turbine1.rho ever became zero, the whole thing would come crashing down from a mathematical point of view. However, as a modeller, I know that the density is always going to be positive, so I don't care about that.

I somehow have the feeling that the implementation in OMC is trying to be too cautious to avoid potential singularities that will never materialize, at the cost of completely useless dynamic state selections that have a huge cost in terms of numerical performance.

Do I miss something here?

That seems to be the case, as i already stated, the static state selection was too strict for StateSelect.never state candidates in the past, therefore it might be the case that we also need to force some other states. IP_turbine1.rho does not seem to have any StateSelect attribute, so i am not quite sure how the compiler should deduce any information about it not reaching zero.

[Karim.Abdelhak]

Replying to casella:

I just wonder: some of the components in this model have (appropriate) StateSelect.prefer attributes, set (if I am not mistaken) on variables that actually appear under derivative operator. Some other components have differentiated variables that should actually be picked as states, but they do not have an explicit StateSelect.prefer attribute.

A fixed choice of states including the ones with StateSelect.prefer and the ones which appear under derivative operator with no StateSelect indication would be just fine, and is in fact what Dymola does. Could it be that the partial indication of StateSelect somehow confuses the algorithm?

It shouldn't confuse it, but maybe we just handle it not as intended. I have an idea on what we could change, but that would take some time to implement and test. And i don't even know if that causes this particular problem.

[Karim.Abdelhak]

Replying to casella:

@Francesco Do you think alias elimination is involved here, or do you think these equations should not be part of the index reduction problem at all?

Most definitely they shouldn't, see above.

From what i can see it may be the case, that we might still analyze the simplified expression at some point during causalization such that it leads to this error.

OMG. All instances of homotopy should be replaced by the actual expression before even starting to tackle causalization. Isn't this the case already?

Homotopy is exclusively relevant for the handling of the initialization problem, which I understand is considered after causalization, when you know who the states are. The backend should restart from scratch at this point, by joining the set of the original equations (including homotopy) with the set of initial equations. Then, a further version of the problem for the lambda = 0 case should be generated, in which all the structural simplifications enabled by that condition are exploited, to facilitate the convergence of the nonlinear solvers at the first step when lambda = 0.

Do I miss something?

It is not actually replaced, but all modules that find a homotopy operator outside of the initialization block should only consider the actual expression. The initialization system is created from the simulation system after Index Reduction, therefore we need to preserve that information for at least that long. Replacing it and replacing it again seems unnecessary, but as i said it might be that not all modules actually do this, hopefully they do but i am not sure.

Some variables like IP_turbine1.wnom only appear in the simplified expression and should not be considered for matching.

IP_turbine1.wnom is a parameter, not a variable. If I am not mistaken, in this case the simplified equations involve a strict subset of the unknowns of the actual problem.

Hope this helps :)

Yes it does :) i think those equations just don't belong in the MSSS therefore i will try to find the reason they got there in the first place.

[casella]

Replying to Karim.Abdelhak:

The equations i posted are only a subset of the full MSSS.

Aha.

The full thing is too large therefore i only post the equation indices to showcase this:

1233, 1175, 1177, 1262, 1171, 56, 29, 55, 30, 36, 31, 32, 48, 2, 49, 1172, 1182, 1191, 1200, 1209, 38, 43, 45, 47, 41, 42, 40, 39, 46, 1180, 1183, 1218, 1219, 1239, 1216, 1217, 1224, 1214, 1215, 1189, 1192, 1245, 1198, 1201, 1251, 1207, 1210, 1257, 1893, 1840, 1866, 1867, 1887, 1864, 1865, 1837, 1872, 1862, 1863, 1839, 1910, 18     19, 3, 65, 51, 50, 1820, 1825, 1830, 1848, 1857, 66, 59, 60, 57, 61, 58, 62, 64, 67, 1823, 1881, 1828, 1831, 1846, 1849, 1899, 1855, 1858, 1905

Wow, that's about 100 variables. No aliases there?

From these equation we try to find a matching to linear appearing states such that they can be chosen statically as dummy states, all remaining variables will be continuous states. The subset i posted cannot be solved for enough linear appearing variables and the compiler will therefore decide to chose the remaining dummy states dynamically (and therefore also the continuous ones dynamically). The rest can be fully chosen statically.

Hmm, this is too terse from me. Are you referring to standard pratice in Sven Erik's dummy derivative algorithm, or is this omc-specific? In the former case, can you refer me to the place in his paper where this is explained?

This static selection has proven to be too strict in the past, therefore i adapted it such that it forces all StateSelect.never states to be matched here, hence to be dummy states, even if they only appear nonlinear.

I'm afraid I still don't grasp the full story, but it seems to me that this is way too strict in general. People very seldom use StateSelect.never, except for debugging or out of sheer desperation. Most of the time it's StateSelect.default or StateSelect.prefer.

Maybe we should come up with some simple conceptual example to understand what we mean.

As a modeller, I am perfectly fine with that. For example, if IP_turbine1.rho ever became zero, the whole thing would come crashing down from a mathematical point of view. However, as a modeller, I know that the density is always going to be positive, so I don't care about that.

I somehow have the feeling that the implementation in OMC is trying to be too cautious to avoid potential singularities that will never materialize, at the cost of completely useless dynamic state selections that have a huge cost in terms of numerical performance.

Do I miss something here?

That seems to be the case, as i already stated, the static state selection was too strict for StateSelect.never state candidates in the past, therefore it might be the case that we also need to force some other states. IP_turbine1.rho does not seem to have any StateSelect attribute, so i am not quite sure how the compiler should deduce any information about it not reaching zero.

Let me put it this way: trying to make sure that the compiler will always find the best state selection automatically by trying to avoid potentially singular nonlinear dependencies at all costs via dynamic state selection is an interesting objective, but I have the feeling that in general it leads to very expensive (and in practice unnecessary) dynamic choice of states.

What I expect as a modeller is that variables with StateSelect.prefer are preferably used as states, as well as variables appearing under derivative sign in the code, unless they are marked with StateSelect.avoid. In general, I don't expect that algebraic variables are selected as states (dynamic or static) unless I explicitly mark them with StateSelect.prefer or StateSelect.always.

What I understand (though I lack some details that we may discuss live) is that the currently implemented criterion is not shy of doing that, based on a structural analysis including linear/nonlinear dependencies, if it thinks it can save the system from some singularity. As far as I understand, that's not what modellers expect from a Modelica tool, though this may be my own interpretation. I always reasoned with these guidelines in my mind:

• If there are only StateSelect.default attributes, the tool will try to use the variables under derivative operator as states, and kick out some of them in case of algebraic constraints between them
• If there are some StateSelect.prefer attributes, the tool should try hard to get them as states. If they are redundant (e.g. in the case of quaternions for free-flying multibody systems) and no specific choice is guaranteed to always be nonsingular, dynamic state selection is applied among them
• If there are someStateSelect.always attributes the tools tries to get them a static states, and if it can't (e.g. if we had StateSelect.always on all four quaternions), it fails
• If there are some StateSelect.avoid or never attributes on variables appearing under derivative operator, the tool should try hard to remove them from the state set. In case of never, it will fail if it can't.
• Under no circumstances should an algebraic variable without StateSelect.prefer/always be promoted to a state, static or dynamic.

I have no idea what is implemented in Dymola (that's a trade secret), but all applications I've dealt with so far are empirically consistent with these guidelines, and I always found the results appropriate. I only got unexpected dynamic state selection when there was some error in the model, and in fact I use that outcome as an indicator of potential issues in the model. This is of course the case if there are no bogus StateSelect.prefer attributes on algebraic variables as in the pump of the HeatingSystem example, which we now removed because it was really misleading.

Modellers are not expecting the tool to outsmart them in the choice of state variables. The reason for this is that doing it properly would require a lot of extra information (e.g., that rho will never become zero) to be explicitly declared. This would be very tedious and it would be hard to come up with complete specifications involving all required properties to avoid an overly strict selection of dummies.

To the contrary, modellers have a good understanding based on their experience and a-priori knowledge of what good choices of states are. IMHO we should rely on this knowledge. By default those are the variables that appear under derivative operators. In case they aren't, they can be marked with StateSelect.avoid/never and/or some algebraic variables can be marked with StateSelect.prefer/always to replace them. Then the tool should try to prefer the preferred and avoid the avoided, while making sure that index 1 is eventually achieved. Modellers never expect the tool to choose algebraic variables other than the ones marked with StateSelect.prefer/always to become states, either statically or dynamically.

Of course I'm not claiming that my view is the only possible one, but that's what I have consistently seen in my experience.

Maybe we should discuss this within the MAP-LANG group. The current specification is extremely vague, but modellers (like me) may have stronger a-priori assumptions with regards to the state selection mechanism, and then it would be good to make them more explicit. Of course this may be dismissed as "tool issue", but it is not; in fact, it is a modelling issue.

Does this list of guidelines help you in getting an algoritm which is less strict and ends up in less dynamic state selections? I guess if you remove a priori from the set of potential states all the algebraic variables that have StateSelect.default/avoid/never everything becomes much simpler, because most combinations you need to deal with now are ruled out a priori.

Of course you could keep the "strict" algorithm we have now as optional, and it may help in some situations where lack of experience can lead to poor state selection. However, I think in most cases an algorithm that follows those guidelines would just be fine.

[casella]

Summarizing today's discussion, @Karim can you try to set the StateSelect attribute of all the variables which are not originally under derivative operator and do not have StateSelect.always/prefer to StateSelect.never, and see what is the impact on the testsuite.

I expect the outcome to be positive, and maybe this could allow to avoid worrying about the correctness of generated Jacobians in these cases.

Thanks!

[casella]

@Karim, sorry to push you on this, but I depend on this model to work for my talk at the Annual Workshop in 1 1/2 months from now. I would be more relaxed if we got this to work before the Christmas break.

This fix seems an easy one, if there are issues with other models in the testsuite but some benefits here you could make it optional with a debug flag for the time being.

Thanks!

[Karim.Abdelhak]

Replying to casella:

@Karim, sorry to push you on this, but I depend on this model to work for my talk at the Annual Workshop in 1 1/2 months from now. I would be more relaxed if we got this to work before the Christmas break.

This fix seems an easy one, if there are issues with other models in the testsuite but some benefits here you could make it optional with a debug flag for the time being.

Thanks!

I tried a lot of things and it all came down to the same set of equations i mentioned in comment:3. I was able to generate a static state selection but the simulation and compilation time seems quite high. The output:

record SimulationResult                
    resultFile = "/home/kab/Repo/FH/promotion/modelicaLibs/ticketTests/5727/Test_res.mat",
    simulationOptions = "startTime = 0.0, stopTime = 1000.0, numberOfIntervals = 500, tolerance = 1e-06, method = '
dassl', fileNamePrefix = 'Test', options = '', outputFormat = 'mat', variableFilter = '.*', cflags = '', simflags =
 ''",               
    messages = "LOG_SUCCESS       | info    | The initialization finished successfully without homotopy method.
LOG_SUCCESS       | info    | The simulation finished successfully.
",                           
    timeFrontend = 1.283448335,                   
    timeBackend = 18.144295305,       
    timeSimCode = 1.385021238,
    timeTemplates = 3.139861082,    
    timeCompile = 31.983834956,
    timeSimulation = 48.052460234,
    timeTotal = 103.989041067                            
end SimulationResult;                                                                                              "Warning: The model contains alias variables with conflicting start and/or nominal values. It is recommended to resolve the conflicts, because otherwise the system could be hard to solve. To print the conflicting alias sets and the chosen candidates please use -d=aliasConflicts.                                                                  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").                        Notification: Model statistics after passing the back-end for initialization:                                       * Number of independent subsystems: 230                                                                            * Number of states: 0 ()                                                                                           * Number of discrete variables: 143 ('-d=discreteinfo' for list of discrete vars)                                  * Number of discrete states: 0 ('-d=discreteinfo' for list of discrete states)
 * Top-level inputs: 0                                                                                             Notification: Strong component statistics for initialization (911):                    
 * Single equations (assignments): 904                                                                              * Array equations: 0                                                          
 * Algorithm blocks: 0                                                                                              * Record equations: 3                                                         
 * When equations: 0                                                                                                * If-equations: 0                                                             
 * Equation systems (linear and non-linear blocks): 0
 * Torn equation systems: 4                                   
 * Mixed (continuous/discrete) equation systems: 0
Notification: Torn system details for strict tearing set:
 * Linear torn systems: 0
 * Non-linear torn systems: 4 {1 16,2 22,1 16,157 1606}
Notification: Model statistics after passing the back-end for simulation:
 * Number of independent subsystems: 7
 * Number of states: 151 (constantSpeed.phi,valveVap1.x,valveVap2.x,valveVap3.x,fWPH_12_1.ValveDynamic.y,fWPH_12_1.
Level_PID.I,fWPH_12_1.flowFW_des.p,fWPH_12_1.flowFW_des.htilde[1],fWPH_12_1.flowFW_des.htilde[2],fWPH_12_1.flowFW_d
es.htilde[3],fWPH_12_1.flowFW_des.htilde[4],fWPH_12_1.flowFW_des.htilde[5],fWPH_12_1.flowFW_cond.p,fWPH_12_1.flowFW
_cond.htilde[1],fWPH_12_1.flowFW_cond.htilde[2],fWPH_12_1.flowFW_cond.htilde[3],fWPH_12_1.flowFW_cond.htilde[4],fWP
H_12_1.flowFW_cond.htilde[5],fWPH_12_1.flowFW_sub.p,fWPH_12_1.flowFW_sub.htilde[1],fWPH_12_1.flowFW_sub.htilde[2],f
WPH_12_1.flowFW_sub.htilde[3],fWPH_12_1.flowFW_sub.htilde[4],fWPH_12_1.flowFW_sub.htilde[5],fWPH_12_1.flowDRAIN.p,f
WPH_12_1.flowDRAIN.htilde[1],fWPH_12_1.flowDRAIN.htilde[2],fWPH_12_1.flowDRAIN.htilde[3],fWPH_12_1.flowDRAIN.htilde
[4],fWPH_12_1.flowDRAIN.htilde[5],fWPH_12_1.flowSTEAM.htilde[1],fWPH_12_1.flowSTEAM.htilde[2],fWPH_12_1.flowSTEAM.h
tilde[3],fWPH_12_1.flowSTEAM.htilde[4],fWPH_12_1.flowSTEAM.htilde[5],fWPH_12_1.metalTubeFV_des.Tvol[1],fWPH_12_1.me
talTubeFV_des.Tvol[2],fWPH_12_1.metalTubeFV_des.Tvol[3],fWPH_12_1.metalTubeFV_des.Tvol[4],fWPH_12_1.metalTubeFV_des
.Tvol[5],fWPH_12_1.metalTubeFV_cond.Tvol[1],fWPH_12_1.metalTubeFV_cond.Tvol[2],fWPH_12_1.metalTubeFV_cond.Tvol[3],f
WPH_12_1.metalTubeFV_cond.Tvol[4],fWPH_12_1.metalTubeFV_cond.Tvol[5],fWPH_12_1.metalTubeFV_sub.Tvol[1],fWPH_12_1.me
talTubeFV_sub.Tvol[2],fWPH_12_1.metalTubeFV_sub.Tvol[3],fWPH_12_1.metalTubeFV_sub.Tvol[4],fWPH_12_1.metalTubeFV_sub
.Tvol[5],fWPH_12_1.nusseltCondenser.zl,fWPH_12_1.nusseltCondenser.hv,fWPH_12_1.nusseltCondenser.hl,fWPH_13_1.ValveD
ynamic.y,fWPH_13_1.Level_PID.I,fWPH_13_1.flowFW_des.p,fWPH_13_1.flowFW_des.htilde[1],fWPH_13_1.flowFW_des.htilde[2]
,fWPH_13_1.flowFW_des.htilde[3],fWPH_13_1.flowFW_des.htilde[4],fWPH_13_1.flowFW_des.htilde[5],fWPH_13_1.flowFW_cond
.p,fWPH_13_1.flowFW_cond.htilde[1],fWPH_13_1.flowFW_cond.htilde[2],fWPH_13_1.flowFW_cond.htilde[3],fWPH_13_1.flowFW
_cond.htilde[4],fWPH_13_1.flowFW_cond.htilde[5],fWPH_13_1.flowFW_sub.p,fWPH_13_1.flowFW_sub.htilde[1],fWPH_13_1.flo
wFW_sub.htilde[2],fWPH_13_1.flowFW_sub.htilde[3],fWPH_13_1.flowFW_sub.htilde[4],fWPH_13_1.flowFW_sub.htilde[5],fWPH
_13_1.flowDRAIN.p,fWPH_13_1.flowDRAIN.htilde[1],fWPH_13_1.flowDRAIN.htilde[2],fWPH_13_1.flowDRAIN.htilde[3],fWPH_13
_1.flowDRAIN.htilde[4],fWPH_13_1.flowDRAIN.htilde[5],fWPH_13_1.flowSTEAM.htilde[1],fWPH_13_1.flowSTEAM.htilde[2],fW
PH_13_1.flowSTEAM.htilde[3],fWPH_13_1.flowSTEAM.htilde[4],fWPH_13_1.flowSTEAM.htilde[5],fWPH_13_1.metalTubeFV_des.T
vol[1],fWPH_13_1.metalTubeFV_des.Tvol[2],fWPH_13_1.metalTubeFV_des.Tvol[3],fWPH_13_1.metalTubeFV_des.Tvol[4],fWPH_1
3_1.metalTubeFV_des.Tvol[5],fWPH_13_1.metalTubeFV_cond.Tvol[1],fWPH_13_1.metalTubeFV_cond.Tvol[2],fWPH_13_1.metalTu
beFV_cond.Tvol[3],fWPH_13_1.metalTubeFV_cond.Tvol[4],fWPH_13_1.metalTubeFV_cond.Tvol[5],fWPH_13_1.metalTubeFV_sub.T
vol[1],fWPH_13_1.metalTubeFV_sub.Tvol[2],fWPH_13_1.metalTubeFV_sub.Tvol[3],fWPH_13_1.metalTubeFV_sub.Tvol[4],fWPH_1
3_1.metalTubeFV_sub.Tvol[5],fWPH_13_1.nusseltCondenser.zl,fWPH_13_1.nusseltCondenser.hv,fWPH_13_1.nusseltCondenser.
hl,fWPH_14_1.ValveDynamic.y,fWPH_14_1.Level_PID.I,fWPH_14_1.flowFW_des.p,fWPH_14_1.flowFW_des.htilde[1],fWPH_14_1.f
lowFW_des.htilde[2],fWPH_14_1.flowFW_des.htilde[3],fWPH_14_1.flowFW_des.htilde[4],fWPH_14_1.flowFW_des.htilde[5],fW
PH_14_1.flowFW_cond.p,fWPH_14_1.flowFW_cond.htilde[1],fWPH_14_1.flowFW_cond.htilde[2],fWPH_14_1.flowFW_cond.htilde[
3],fWPH_14_1.flowFW_cond.htilde[4],fWPH_14_1.flowFW_cond.htilde[5],fWPH_14_1.flowFW_sub.p,fWPH_14_1.flowFW_sub.htil
de[1],fWPH_14_1.flowFW_sub.htilde[2],fWPH_14_1.flowFW_sub.htilde[3],fWPH_14_1.flowFW_sub.htilde[4],fWPH_14_1.flowFW
_sub.htilde[5],fWPH_14_1.flowDRAIN.p,fWPH_14_1.flowDRAIN.htilde[1],fWPH_14_1.flowDRAIN.htilde[2],fWPH_14_1.flowDRAI
N.htilde[3],fWPH_14_1.flowDRAIN.htilde[4],fWPH_14_1.flowDRAIN.htilde[5],fWPH_14_1.flowSTEAM.htilde[1],fWPH_14_1.flo
wSTEAM.htilde[2],fWPH_14_1.flowSTEAM.htilde[3],fWPH_14_1.flowSTEAM.htilde[4],fWPH_14_1.flowSTEAM.htilde[5],fWPH_14_
1.metalTubeFV_des.Tvol[1],fWPH_14_1.metalTubeFV_des.Tvol[2],fWPH_14_1.metalTubeFV_des.Tvol[3],fWPH_14_1.metalTubeFV
_des.Tvol[4],fWPH_14_1.metalTubeFV_des.Tvol[5],fWPH_14_1.metalTubeFV_cond.Tvol[1],fWPH_14_1.metalTubeFV_cond.Tvol[2
],fWPH_14_1.metalTubeFV_cond.Tvol[3],fWPH_14_1.metalTubeFV_cond.Tvol[4],fWPH_14_1.metalTubeFV_cond.Tvol[5],fWPH_14_
1.metalTubeFV_sub.Tvol[1],fWPH_14_1.metalTubeFV_sub.Tvol[2],fWPH_14_1.metalTubeFV_sub.Tvol[3],fWPH_14_1.metalTubeFV
_sub.Tvol[4],fWPH_14_1.metalTubeFV_sub.Tvol[5],fWPH_14_1.nusseltCondenser.zl,fWPH_14_1.nusseltCondenser.hv,fWPH_14_
1.nusseltCondenser.hl)
 * Number of discrete variables: 292 ('-d=discreteinfo' for list of discrete vars)
 * Number of discrete states: 0 ('-d=discreteinfo' for list of discrete states)
 * Top-level inputs: 0
Notification: Strong component statistics for simulation (1781):
 * Single equations (assignments): 1515
 * Array equations: 0
 * Algorithm blocks: 0
 * Record equations: 138
 * When equations: 0
 * If-equations: 0
 * Equation systems (linear and non-linear blocks): 1
 * Torn equation systems: 127
 * Mixed (continuous/discrete) equation systems: 0
Notification: Equation system details:
 * Constant Jacobian: 0
 * Linear Jacobian (size,density): 0
 * Non-linear Jacobian: 1 {1}
 * Without analytic Jacobian: 0
Notification: Torn system details for strict tearing set:
 * Linear torn systems: 87 {(6,72.2%) 10,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(6,72.2%)
 10,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(6,72.2%) 10,(1,100.0%) 1,(1,100.0%) 1,(1,100.
0%) 1,(1,100.0%) 1,(1,100.0%) 1,(6,72.2%) 10,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(6,72
.2%) 10,(1,100.0%) 2,(6,72.2%) 10,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(6,72.2%) 10,(1,
100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(
6,72.2%) 10,(1,100.0%) 2,(6,72.2%) 10,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2
,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(6,72.2%) 10,(1,100.0%) 2,(6,72.2%) 10,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%)
 2,(1,100.0%) 2,(1,100.0%) 2,(6,72.2%) 10,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0
%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100
.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,1
00.0%) 1}
 * Non-linear torn systems: 40 {1 15,2 24,9 20,1 5,9 20,1 5,1 15,9 20,1 5,1 8,1 8,1 8,1 8,1 8,1 8,1 8,1 8,1 1,1 1,1
 1,1 1,1 1,1 33,1 5,1 1,1 1,1 1,1 1,1 1,1 33,1 5,1 1,1 1,1 1,1 1,1 1,1 33,3 18,1 2,1 2}
[/home/kab/Repo/FH/promotion/modelicaLibs/ticketTests/5727/Test.mo:1315:7-1315:35:writable] Warning: Variable valve
Vap3.x has attribute stateSelect=StateSelect.never, but was selected as a state
[/home/kab/Repo/FH/promotion/modelicaLibs/ticketTests/5727/Test.mo:1315:7-1315:35:writable] Warning: Variable valve
Vap2.x has attribute stateSelect=StateSelect.never, but was selected as a state
[/home/kab/Repo/FH/promotion/modelicaLibs/ticketTests/5727/Test.mo:1315:7-1315:35:writable] Warning: Variable valve
Vap1.x has attribute stateSelect=StateSelect.never, but was selected as a state

As you can see some of the StateSelect.never vars could not be forced in, i guess these are not user tagged but rather algebraic states that got into state selection.

I think i can do a little something about the compilation time since right now i am full matching the never vars to get perfect information, but i can reuse old information here. Won't be much of a save but still a little.

Also i need to check if it is viable with the testsuite. Is that about what you expected @Francesco?

[casella]

Replying to Karim.Abdelhak:

Also i need to check if it is viable with the testsuite. Is that about what you expected @Francesco?

Unfortunately not at all. Dymola performs a static state selection and simulates quickly and without any problem. We are probably overlooking something.

[casella]

@Karim

I went through the list of the 151 statically selected states reported in comment:13 and compared it with Dymola's (which is also what I would expect).

The list is the same, except for Dymola selecting some aliases (h[2:6] instead of htilde[1:5]), which is irrelevant, and except these three variables:

Dymola:

• fWPH_12_1.nusseltCondenser.p
• fWPH_13_1.nusseltCondenser.p
• fWPH_13_1.nusseltCondenser.p

OMC:

• valveVap1.x
• valveVap2.x
• valveVap3.x

This is really weird, because NusselCondenser.p has StateSelect.prefer, while ValveVap.x has StateSelect.default; besides, ValveVap.x does not appear under derivative sign, so it should actually get StateSelect.never according to my proposal in comment:11

Can you pinpoint why your algorithm is selecting the latter instead of the former? This really makes no sense to me.

[casella]

@Karim, maybe I've got an important clue for you.

The ValveVap model contains the equation

x = dp/p;

where p is the inlet pressure of the valve, which is equal to NusseltCondenser.p by way of a connect statement.

It looks like your algorithm is scared about this potential division by zero, so if you really force a static state selection it prefers using ValveVap.x instead of NusseltCondenser.p, in order to obtain a safer one, nonwithstanding the indications from the StateSelect attributes.

This is where my argument of comment:10 enters into play. As a modeller, I know that the condenser pressure will never become zero. In fact, the model will fail much earlier than that, because of asserts on lower pressure limits on the IF97 steam model. So, it is perfectly safe to use it as a statically chosen state, even if it appears at the denominator in some equation. People have used absolute pressures as states in models for decades and nobody ever got hurt by that :)

In other words, user-defined StateSelect.prefer attributes should override the negative outcome of any checks about nonlinear dependencies stemming from that choice o states. Only in the case there are redundant StateSelect.prefer potential states one should check that no singularity occurs (e.g. as in the case of free-flying bodies), but one should not rely on redundancy given by algebraic variables without an explicitly declared StateSelect.prefer/always attribute.

Hope this helps.

[Karim.Abdelhak]

I did some changes and ended up with the changes you expected. Basically i changed a little bit in the matching such that it does not greedily choose but instead respects the presorted order according to StateSelect attributes. Unfortunately we still have an overall time of 100 secs.

The new state choice:

resultFile = "/home/kab/Repo/FH/promotion/modelicaLibs/ticketTests/5727/Test_res.mat",               [34/92508]
    simulationOptions = "startTime = 0.0, stopTime = 1000.0, numberOfIntervals = 500, tolerance = 1e-06, method = '
dassl', fileNamePrefix = 'Test', options = '', outputFormat = 'mat', variableFilter = '.*', cflags = '', simflags =
 ''",               
    messages = "LOG_SUCCESS       | info    | The initialization finished successfully without homotopy method.
LOG_SUCCESS       | info    | The simulation finished successfully.
",                           
    timeFrontend = 1.407411655,                   
    timeBackend = 16.119984532,       
    timeSimCode = 1.385240399,
    timeTemplates = 3.19501146,     
    timeCompile = 31.935285696,
    timeSimulation = 45.336503616,
    timeTotal = 99.37955881400001                        
end SimulationResult;                                                                                              "Warning: The model contains alias variables with conflicting start and/or nominal values. It is recommended to resolve the conflicts, because otherwise the system could be hard to solve. To print the conflicting alias sets and the chosen candidates please use -d=aliasConflicts.                                                                  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").                        Notification: Model statistics after passing the back-end for initialization:                                       * Number of independent subsystems: 230                                                                            * Number of states: 0 ()                                                                                           * Number of discrete variables: 143 ('-d=discreteinfo' for list of discrete vars)                                  * Number of discrete states: 0 ('-d=discreteinfo' for list of discrete states)
 * Top-level inputs: 0                                                                                             Notification: Strong component statistics for initialization (911):
 * Single equations (assignments): 904
 * Array equations: 0                                                            
 * Algorithm blocks: 0
 * Record equations: 3
 * When equations: 0
 * If-equations: 0
 * Equation systems (linear and non-linear blocks): 0
 * Torn equation systems: 4                                                                                   
 * Mixed (continuous/discrete) equation systems: 0            
Notification: Torn system details for strict tearing set:
 * Linear torn systems: 0
 * Non-linear torn systems: 4 {1 16,2 22,1 16,157 1606}
Notification: Model statistics after passing the back-end for simulation:
 * Number of independent subsystems: 10
 * Number of states: 151 (constantSpeed.phi,fWPH_12_1.ValveDynamic.y,fWPH_12_1.Level_PID.I,fWPH_12_1.flowFW_des.p,f
WPH_12_1.flowFW_des.htilde[1],fWPH_12_1.flowFW_des.htilde[2],fWPH_12_1.flowFW_des.htilde[3],fWPH_12_1.flowFW_des.ht
ilde[4],fWPH_12_1.flowFW_des.htilde[5],fWPH_12_1.flowFW_cond.p,fWPH_12_1.flowFW_cond.htilde[1],fWPH_12_1.flowFW_con
d.htilde[2],fWPH_12_1.flowFW_cond.htilde[3],fWPH_12_1.flowFW_cond.htilde[4],fWPH_12_1.flowFW_cond.htilde[5],fWPH_12
_1.flowFW_sub.p,fWPH_12_1.flowFW_sub.htilde[1],fWPH_12_1.flowFW_sub.htilde[2],fWPH_12_1.flowFW_sub.htilde[3],fWPH_1
2_1.flowFW_sub.htilde[4],fWPH_12_1.flowFW_sub.htilde[5],fWPH_12_1.flowDRAIN.p,fWPH_12_1.flowDRAIN.htilde[1],fWPH_12
_1.flowDRAIN.htilde[2],fWPH_12_1.flowDRAIN.htilde[3],fWPH_12_1.flowDRAIN.htilde[4],fWPH_12_1.flowDRAIN.htilde[5],fW
PH_12_1.flowSTEAM.htilde[1],fWPH_12_1.flowSTEAM.htilde[2],fWPH_12_1.flowSTEAM.htilde[3],fWPH_12_1.flowSTEAM.htilde[
4],fWPH_12_1.flowSTEAM.htilde[5],fWPH_12_1.metalTubeFV_des.Tvol[1],fWPH_12_1.metalTubeFV_des.Tvol[2],fWPH_12_1.meta
lTubeFV_des.Tvol[3],fWPH_12_1.metalTubeFV_des.Tvol[4],fWPH_12_1.metalTubeFV_des.Tvol[5],fWPH_12_1.metalTubeFV_cond.
Tvol[1],fWPH_12_1.metalTubeFV_cond.Tvol[2],fWPH_12_1.metalTubeFV_cond.Tvol[3],fWPH_12_1.metalTubeFV_cond.Tvol[4],fW
PH_12_1.metalTubeFV_cond.Tvol[5],fWPH_12_1.metalTubeFV_sub.Tvol[1],fWPH_12_1.metalTubeFV_sub.Tvol[2],fWPH_12_1.meta
lTubeFV_sub.Tvol[3],fWPH_12_1.metalTubeFV_sub.Tvol[4],fWPH_12_1.metalTubeFV_sub.Tvol[5],fWPH_12_1.nusseltCondenser.
zl,fWPH_12_1.nusseltCondenser.p,fWPH_12_1.nusseltCondenser.hv,fWPH_12_1.nusseltCondenser.hl,fWPH_13_1.ValveDynamic.
y,fWPH_13_1.Level_PID.I,fWPH_13_1.flowFW_des.p,fWPH_13_1.flowFW_des.htilde[1],fWPH_13_1.flowFW_des.htilde[2],fWPH_1
3_1.flowFW_des.htilde[3],fWPH_13_1.flowFW_des.htilde[4],fWPH_13_1.flowFW_des.htilde[5],fWPH_13_1.flowFW_cond.p,fWPH
_13_1.flowFW_cond.htilde[1],fWPH_13_1.flowFW_cond.htilde[2],fWPH_13_1.flowFW_cond.htilde[3],fWPH_13_1.flowFW_cond.h
tilde[4],fWPH_13_1.flowFW_cond.htilde[5],fWPH_13_1.flowFW_sub.p,fWPH_13_1.flowFW_sub.htilde[1],fWPH_13_1.flowFW_sub
.htilde[2],fWPH_13_1.flowFW_sub.htilde[3],fWPH_13_1.flowFW_sub.htilde[4],fWPH_13_1.flowFW_sub.htilde[5],fWPH_13_1.f
lowDRAIN.p,fWPH_13_1.flowDRAIN.htilde[1],fWPH_13_1.flowDRAIN.htilde[2],fWPH_13_1.flowDRAIN.htilde[3],fWPH_13_1.flow
DRAIN.htilde[4],fWPH_13_1.flowDRAIN.htilde[5],fWPH_13_1.flowSTEAM.htilde[1],fWPH_13_1.flowSTEAM.htilde[2],fWPH_13_1
.flowSTEAM.htilde[3],fWPH_13_1.flowSTEAM.htilde[4],fWPH_13_1.flowSTEAM.htilde[5],fWPH_13_1.metalTubeFV_des.Tvol[1],
fWPH_13_1.metalTubeFV_des.Tvol[2],fWPH_13_1.metalTubeFV_des.Tvol[3],fWPH_13_1.metalTubeFV_des.Tvol[4],fWPH_13_1.met
alTubeFV_des.Tvol[5],fWPH_13_1.metalTubeFV_cond.Tvol[1],fWPH_13_1.metalTubeFV_cond.Tvol[2],fWPH_13_1.metalTubeFV_co
nd.Tvol[3],fWPH_13_1.metalTubeFV_cond.Tvol[4],fWPH_13_1.metalTubeFV_cond.Tvol[5],fWPH_13_1.metalTubeFV_sub.Tvol[1],
fWPH_13_1.metalTubeFV_sub.Tvol[2],fWPH_13_1.metalTubeFV_sub.Tvol[3],fWPH_13_1.metalTubeFV_sub.Tvol[4],fWPH_13_1.met
alTubeFV_sub.Tvol[5],fWPH_13_1.nusseltCondenser.zl,fWPH_13_1.nusseltCondenser.p,fWPH_13_1.nusseltCondenser.hv,fWPH_
13_1.nusseltCondenser.hl,fWPH_14_1.ValveDynamic.y,fWPH_14_1.Level_PID.I,fWPH_14_1.flowFW_des.p,fWPH_14_1.flowFW_des
.htilde[1],fWPH_14_1.flowFW_des.htilde[2],fWPH_14_1.flowFW_des.htilde[3],fWPH_14_1.flowFW_des.htilde[4],fWPH_14_1.f
lowFW_des.htilde[5],fWPH_14_1.flowFW_cond.p,fWPH_14_1.flowFW_cond.htilde[1],fWPH_14_1.flowFW_cond.htilde[2],fWPH_14
_1.flowFW_cond.htilde[3],fWPH_14_1.flowFW_cond.htilde[4],fWPH_14_1.flowFW_cond.htilde[5],fWPH_14_1.flowFW_sub.p,fWP
H_14_1.flowFW_sub.htilde[1],fWPH_14_1.flowFW_sub.htilde[2],fWPH_14_1.flowFW_sub.htilde[3],fWPH_14_1.flowFW_sub.htil
de[4],fWPH_14_1.flowFW_sub.htilde[5],fWPH_14_1.flowDRAIN.p,fWPH_14_1.flowDRAIN.htilde[1],fWPH_14_1.flowDRAIN.htilde
[2],fWPH_14_1.flowDRAIN.htilde[3],fWPH_14_1.flowDRAIN.htilde[4],fWPH_14_1.flowDRAIN.htilde[5],fWPH_14_1.flowSTEAM.h
tilde[1],fWPH_14_1.flowSTEAM.htilde[2],fWPH_14_1.flowSTEAM.htilde[3],fWPH_14_1.flowSTEAM.htilde[4],fWPH_14_1.flowST
EAM.htilde[5],fWPH_14_1.metalTubeFV_des.Tvol[1],fWPH_14_1.metalTubeFV_des.Tvol[2],fWPH_14_1.metalTubeFV_des.Tvol[3]
,fWPH_14_1.metalTubeFV_des.Tvol[4],fWPH_14_1.metalTubeFV_des.Tvol[5],fWPH_14_1.metalTubeFV_cond.Tvol[1],fWPH_14_1.m
etalTubeFV_cond.Tvol[2],fWPH_14_1.metalTubeFV_cond.Tvol[3],fWPH_14_1.metalTubeFV_cond.Tvol[4],fWPH_14_1.metalTubeFV
_cond.Tvol[5],fWPH_14_1.metalTubeFV_sub.Tvol[1],fWPH_14_1.metalTubeFV_sub.Tvol[2],fWPH_14_1.metalTubeFV_sub.Tvol[3]
,fWPH_14_1.metalTubeFV_sub.Tvol[4],fWPH_14_1.metalTubeFV_sub.Tvol[5],fWPH_14_1.nusseltCondenser.zl,fWPH_14_1.nussel
tCondenser.p,fWPH_14_1.nusseltCondenser.hv,fWPH_14_1.nusseltCondenser.hl)
 * Number of discrete variables: 292 ('-d=discreteinfo' for list of discrete vars)
 * Number of discrete states: 0 ('-d=discreteinfo' for list of discrete states)
 * Top-level inputs: 0
Notification: Strong component statistics for simulation (1843):
 * Single equations (assignments): 1558
 * Array equations: 0
 * Algorithm blocks: 0
 * Record equations: 162
 * When equations: 0
 * If-equations: 0
 * Equation systems (linear and non-linear blocks): 1
 * Torn equation systems: 122
 * Mixed (continuous/discrete) equation systems: 0
Notification: Equation system details:
 * Constant Jacobian: 0
 * Linear Jacobian (size,density): 0
 * Non-linear Jacobian: 1 {1}
 * Without analytic Jacobian: 0
Notification: Torn system details for strict tearing set:
 * Linear torn systems: 87 {(6,72.2%) 10,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(6,72.2%)
 10,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(6,72.2%) 10,(1,100.0%) 1,(1,100.0%) 1,(1,100.
0%) 1,(1,100.0%) 1,(1,100.0%) 1,(6,72.2%) 10,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(6,72
.2%) 10,(1,100.0%) 2,(6,72.2%) 10,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(6,72.2%) 10,(1,
100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(6,72.2%) 10,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(
1,100.0%) 2,(1,100.0%) 2,(6,72.2%) 10,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(6,72.2%) 10
,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(6,72.2%) 10,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%)
 2,(1,100.0%) 2,(1,100.0%) 2,(6,72.2%) 10,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0
%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 2,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100
.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,100.0%) 1,(1,1
00.0%) 1}
 * Non-linear torn systems: 35 {1 15,2 24,9 21,1 5,9 21,1 5,1 15,9 21,1 5,1 8,1 8,1 8,1 1,1 1,1 1,1 1,1 1,1 8,1 8,1
 8,1 1,1 1,1 1,1 1,1 1,1 8,1 8,1 1,1 1,1 1,1 1,1 1,2 18,2 18,2 29}
"

Let me know if you find anything else that is odd!

[casella]

The choice of states is now ok. If you push your changes I can build a Windows nightly, run the debugger with performance analysis, and try to figure out where all the time is spent.

[casella]

Hint: Dymola reports:

Sizes of nonlinear systems of equations: 
{25, 2, 2, 2, 2, 2, 6, 6, 6, 6, 18, 2, 2, 2, 2, 2, 6, 18, 2, 2, 2, 2, 2, 6, 6, 6, 6, 6, 6, 1}
Sizes after manipulation of the nonlinear systems: 
{2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 2, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1}

OpenModelica reports:

Non-linear torn systems: 35 
{1 15,2 24,9 21,1 5,9 21,1 5,1 15,9 21,1 5,1 8,1 8,1 8,1 1,1 1,1 1,1 1,1 1,1 8,1 8,1 8,1 1,1 1,1 1,1 1,1 1,1 8,1 8,1 1,1 1,1 1,1 1,1 1,2 18,2 18,2 29}

If I interpret this output correcly (each pair of numbers indicates the tearing and torn variables, respectively), we are left with more nonlinear systems of larger size, and after tearing we get three systems with 9 tearing variables, while Dymola gets at most 2 tering variables for each system. I guess this may explain the huge overhead, since I suspect the IF97 water model to be involved in those loops.

Maybe this is because of some missing optimization of IF97 computation (like LateInline). I can only figure this out if I can run the simulation myself with the OMEdit nightly.

[Karim.Abdelhak]

​PR contains the StateSelect.never update,

[casella]

I just got the ​regression report. 11 improvements and 81 regressions.

Big negative impact on the Chemical library. Some models, e.g. ​EnzymeKinetics correctly select dynamic state sets (Dymola also does it), but then fail with

Error: Internal error IndexReduction.pantelidesIndexReduction failed!
System is structurally singulare and cannot handled because number of unassigned
equations is larger than number of states. Use -d=bltdump to get more information.
Error: Internal error Transformation Module PFPlusExt index Reduction Method 
Pantelides failed!

while others, e.g. ​HeatingOfAlcohol fail at initialization. Maybe the wrong initial choice of dynamic states is made, leading to a singularity.

If you could check what changed between this version and the previous one, that would be great.

I'll try the nighly tomorrow on my example to continue the analysis.

[Karim.Abdelhak]

I fixed one small thing which should have a great effect on ModelicaTest.Fluid.TestUtilities.TestRegRoot2Derivatives and similar models. It was just a small error in the update module for differentiation counting.

Furthermore i analyzed ModelicaTest.Fluid.TestExamplesVariants.BranchingDynamicPipes_MomentumSteadyState which has a failure at initialization just as you described. I looked into the generated code and stumbled upon the functions saturationPressure, saturationPressureLiquid and their derivatives in the Media library. It seems like the first function uses the latter one and that is also true for the derivatives. It seems like this inheritance produces problems since the generated function for the derivative of saturationPressureLiquid in c-code seems like utter garbage. It could be some form of automated differentiation result.

I think this is common usage of predefined derivatives of functions, so i don't really understand why it breaks in this particular case. I will have a look at it.

[Karim.Abdelhak]

Some models are fixed now: ​report

[casella]

Replying to Karim.Abdelhak:

Some models are fixed now: ​report

Very good!

[casella]

Update on the original test case. Now that the static selection issue has been fixed, as well as some other issues with the backend, the situation has much improved. Now the simulation time has gone down from 30 s to 3.2 seconds. This is not yet as good as Dymola, which takes about half of that time, but is now perfectly acceptable.

The initialization time with homotopy is still quite long, about 21 s compared to 3.5 seconds, we should also investigate that, but this is out of scope for this ticket

[Karim.Abdelhak]

Replying to casella:

Update on the original test case. Now that the static selection issue has been fixed, as well as some other issues with the backend, the situation has much improved. Now the simulation time has gone down from 30 s to 3.2 seconds. This is not yet as good as Dymola, which takes about half of that time, but is now perfectly acceptable.

The initialization time with homotopy is still quite long, about 21 s compared to 3.5 seconds, we should also investigate that, but this is out of scope for this ticket

Unfortunately i am not even able to reproduce that good behaviour. I still got ~100 seconds total (~60 simulation). :(

[casella]

Given the performance I obtained on this model and similar ones, I think the specific issues addressed by this ticket can be considered as resolved. We still have slow initialization, but that's because of #3921, and we still have slightly sub-optimal simulation, probably due to residual issues with CSE.