Issues with redundant initial conditions in Buildings model?

[Francesco Casella]

Please check ​Buildings.Experimental.DistrictHeatingCooling.Validation.HeatingCoolingHotWaterSmall. The backend fails with

Warning: It was not possible to determine if the initialization problem is consistent, because of not evaluable parameters/start values during compile time: pla.coo.con.vol.dynBal.medium.p = pla.coo.con.vol.dynBal.p_start (pla.sinCoo.p = pla.coo.con.vol.dynBal.p_start)
Warning: It was not possible to determine if the initialization problem is consistent, because of not evaluable parameters/start values during compile time: pla.hea.eva.vol.dynBal.medium.p = pla.hea.eva.vol.dynBal.p_start (pla.sinHea.p = pla.hea.eva.vol.dynBal.p_start)
Notification: Performance of analyzeInitialSystem (initialization): time 0.195/7.507, allocations: 29.7 MB / 2.141 GB, free: 0.611 GB / 1.292 GB
Notification: Performance of solveInitialSystemEqSystem (initialization): time 9.535e-05/7.507, allocations: 24 kB / 2.141 GB, free: 0.611 GB / 1.292 GB
Error: Internal error IndexReduction.pantelidesIndexReduction failed! System is structurally singular and cannot be handled because the number of unassigned equations is larger than the number of states. Use -d=bltdump to get more information.
Error: No system for the symbolic initialization was generated

I guess the problem is that there initial equations made redundant by state constraints, and the backend can't deal with them because they are parameter-dependent.

In this case, I would suggest the backend to make those parameters structural and constant-evaluate them. I tried to do that by turning on the evaluate all parameters option, but unfortunately this doesn't work because the model uses external functions and the frontend cannot evaluate them.

[Andreas Heuermann]

The relevant function getConsistentEquation has a matchcontinue. I'm not a big fan and would replace it directly. The newBackend won't use matchcontinue at all. It makes finding errors very difficult.

With a bit of debugging it should be simple to find the exact location of the thrown error.

Either fixing getConsistentEquation or move the check to simulation time for cases where it is not possible to resolve all parameters at compile time.

[Francesco Casella]

Replying to AnHeuermann:

move the check to simulation time for cases where it is not possible to resolve all parameters at compile time.

We may assume they are all equal for the equation processing (so that getConsistentEquation() can eliminate the redundant ones), and then check this condition at runtime, throwing an error if it doesn't hold. If this turns out to be complicated, just evaluate them at compile time.

Being able to change these parameters without recompiling is of course a nice-to-have feature but it's much less important than getting the model to actually run.

[Francesco Casella]

Ticket retargeted after milestone closed

[Francesco Casella]

Smaller example case with the same issue: ​Buildings.Fluid.HeatExchangers.DXCoils.AirCooled.Examples.MultiStage

[Francesco Casella]

@Karim, maybe I got a good lead on this issue. Please check this example case: ​Buildings.Fluid.HeatExchangers.Validation.HeaterCooler_u.

Diagnostics with -d=bltdump:

[6] 22:35:32 Symbolic Error
[Buildings.Fluid.Interfaces.ConservationEquation: 319:7-319:37]: 
Model is structurally singular, error found sorting equations
  28: $DER.heaDyn.vol.dynBal.medium.p = 0.0
  6: $DER.heaDyn.vol.dynBal.medium.p = 0.0
for variables
  109: $DER.heaDyn.vol.dynBal.medium.p:DUMMY_DER(flow=false fixed = false )  "Absolute pressure of medium"Buildings.Fluid.HeatExchangers.Validation.HeaterCooler_u, Buildings.Fluid.HeatExchangers.HeaterCooler_u, Buildings.Fluid.MixingVolumes.MixingVolume, Buildings.Fluid.Interfaces.ConservationEquation, heaDyn.vol.dynBal.Medium.BaseProperties type: Real
  36: heaDyn.vol.dynBal.ports_mXi_flow[1,1]:VARIABLE(unit = "kg/s" protected = true final = true ) Buildings.Fluid.HeatExchangers.Validation.HeaterCooler_u, Buildings.Fluid.HeatExchangers.HeaterCooler_u, Buildings.Fluid.MixingVolumes.MixingVolume, Buildings.Fluid.Interfaces.ConservationEquation type: Real [2,1]

Now, there are two equations stating that the derivative of a certain pressure is zero, one belonging to the "dynamic" set and the other belonging to the "initial" set.

The initial one comes straight from the initial equations set. I tried to figure out where the dynamic one comes from. If I turn on -d=optdaedump, that equation first appears in the log immediately after "Index reduction done.". It shows up in the previous list, marked "pantelidesIndexReduction (unspecified partition)", as der(heaDyn.vol.dynBal.medium.p) = 0.0. As I understand, this equation comes from the fact that the medium pressure is equal to the fixed sink pressure because of the direct connection between the two volumes, so the Pantelides algorithm identifies this constraint and differentiates it.

Now, the algorithm that identifies redundant initial conditions due to index reduction should detect this pair as such and remove one. For some reason, however, this pair escapes the detection. Can you try to figure out why?

[Francesco Casella]

I tried to build an MWE that recaps the structure of HeaterCooler_u as far as index reduction is concerned:

model M3
  Real T(stateSelect = StateSelect.prefer);
  Real p(stateSelect = StateSelect.prefer);
  Real rho;
  Real u(stateSelect = StateSelect.avoid);
  Real U(stateSelect = StateSelect.avoid);
  Real h;
  Real m(stateSelect = StateSelect.avoid);
  Real w1 = 1 + sin(time);
  Real Q = sin(time);
  Real w2;
  constant Real R = 1;
  constant Real cp = 1;
  constant Real cv = 0.7;
  parameter Real V = 1;
  parameter Real h1 = 1;
  parameter Real p2 = 10;
equation
  m = rho*V;
  rho = p/R*T;
  u = cv*T;
  U = m*u;
  h = cp*T;
  der(m) = w1 + w2;
  der(U) = w1*h1 + w2*h + Q;
  p = p2;
initial equation
  der(p) = 0;
  der(T) = 0;
end M3;

however, in this case the redundant initial equation detection algorithm works fine:

[2] 23:52:13 Translation Warning
The initial conditions are over specified. The following 1 initial equations are redundant, so they are removed from the initialization sytem:
         $DER.p = 0.0.

Why doesn't it for HeaterCooler_u? Any clue?

[Karim Adbdelhak]

I really could not figure out why the current algorithm can't solve this so i implemented an ​entirely new approach. It is able to simulate Buildings.Fluid.HeatExchangers.Validation.HeaterCooler_u, i have not tested the other ones.

We can expect a lot of failures in the testsuite since it is a major change and i just did minor effort to get it running. Let's wait and see what impact it has. I will try to improve it such that it can replace the old stuff, but as i mentioned in our meeting - I have no high hopes.

[Francesco Casella]

OK, good that you found an algorithm to deal with HeaterCooler_u, good luck with the testsuite :)

BTW, I don't see any problem at making multiple attempts with different algorithms, there is no need to have a single jack-of-all-trades algorithm!

What I mean is: if index-reduction takes place (not just state variable change), it is to be expected that we may have redundant initial conditions, unless the modeller is aware of the high-index nature of the problem and manually set the model parameters to get rid of some initial conditions - I consider this as a last resort, but it is not the preferred approach.

So, if the initialization problem was square before index reduction, and becomes overdetermined after it, then we can launch multiple attempts to identify redundant but consistent equations. If one fails, we try another. In fact, one could find a subset of redundant equations and remove them, but not be able to complete the task, so it could pass whatever it obtained to the next one.

Machiavelli's approach: the ends justify the means :)

What do you think?

[Francesco Casella]

I also have another idea, which some friends in Lund may consider heretic, but which could turn out to be extremely effective.

In 99% of the cases I am aware of (mostly thermo-fluid models), the redundant equations of the initialization problem have a very simple form, either der(x) = 0 or x = some_value, where x are state or dummy state variables, and they simply show up more than one time in the final equation set.

Assuming the V-nodes are ordered lexicographically, we just iterate over them, and for each of the ones that correspond to states or derivatives, we check if there is more than one arc connecting to an equation node that has no other arcs. In that case, it is very likely that those equation nodes correspond to redundant initial equations, and we can symbolically check if they are consistent or not.

We could then keep all the other more elaborate algorithms to take care of other potentially occurring cases I'm now aware of.

@Karim , what do you think?

[Karim Adbdelhak]

There are some jenkins issues going on, currently no PR can be tested. I will resume work once it is fixed and we can test pull requests again. ​Current PR: 7903

[Karim Adbdelhak]

Very good news. The tests went really well, only 22 failing tests of which 21 are only due to small equation removal differences(sometimes the redundant equation to remove is arbitrary) or missing consistency check, which is not a big step to implement (probably i will just reuse the old implementation for that).

The only model that actually fails is: /simulation/modelica/msl22/Test3PhaseInitOver.mos. This one is a little intricate since it also reports initial equations which cannot be checked for consistency.

For now i will revive the consistency check and see if it does some magic on the actually failing model.

[Karim Adbdelhak]

There also is one open thing to discuss. The new implementation trys as hard as it can to only fix variables that are quote on quote fixable. (states, discretes, parameters, ..) but if that does not work it also finds the least amount of necessary remaining variables to fix.

For well posed models this will never occur, but it would be able to tryhard on ill posed initialization with that. So my natural questions is: What happens if we have no other choice but fix unfixable variables?

• Option 1: Fail - same behavior as old implementation
• Option 2: Fail but offer a compiler flag -ignoreIllPosedInitialization or something less convoluted
• Option 3: Issue a warning for the illegally fixed variables and continue as if nothing happened.

Keep in mind that there is no extra work in offering that option, the compilation process is just as fast as without it and well behaved models are not affected at all. Should we protect the user from bad decisions here or is it viable to possibly initialize in a bad way rather than not at all?

Furthermore i want to notice that there is no maxMixedIndex anymore with this. Everything is solved in one swoop no matter the index.

[Francesco Casella]

Replying to Karim.Abdelhak:

Very good news. The tests went really well, only 22 failing tests of which 21 are only due to small equation removal differences(sometimes the redundant equation to remove is arbitrary) or missing consistency check, which is not a big step to implement (probably i will just reuse the old implementation for that).

The only model that actually fails is: /simulation/modelica/msl22/Test3PhaseInitOver.mos. This one is a little intricate since it also reports initial equations which cannot be checked for consistency.

For now i will revive the consistency check and see if it does some magic on the actually failing model.

Excellent! Pls keep me posted.

[Francesco Casella]

Replying to Karim.Abdelhak:

There also is one open thing to discuss. The new implementation trys as hard as it can to only fix variables that are quote on quote fixable. (states, discretes, parameters, ..) but if that does not work it also finds the least amount of necessary remaining variables to fix.

For well posed models this will never occur, but it would be able to tryhard on ill posed initialization with that.

@Karim, I'm not sure I get your point. Additional variables need to be fixed if there are not enough initial equations to get a square, well-posed initialization problem. What we are trying to solve here is the opposite problem: what do we do when we have too many equations (because of initial equations and index reduction)?

Are you referring to the case where a part of the model doesn't have enought initial equations, and another one has one too many because of index reduction, so you get a square mix-determined system?

[Karim Adbdelhak]

During my work on the new initialization i activated ASSC just to see if it fixes some problems. We talked about this shortly in our dev-meetings and you asked for an example where this would be even necessary. I stumbled upon this old reoccurring problematic model: Modelica.Mechanics.MultiBody.Examples.Constraints.PrismaticConstraint.mos (currently testing with MSL 3.2 version).

In our old implementation the compiler warns about a singular system during initialization:

Warning: The linear system: 
1 : bodyOfConstraint.v_0[1] = 0.2148655148851041 * freeMotionScalarInit.v_rel_a_1 + (-0.5318108111446446) * freeMotionScalarInit.v_rel_a_2 + 0.8191520442889918 * freeMotionScalarInit.v_rel_a_3
2 : 0.0 = 0.8191520442889918 * bodyOfConstraint.v_0[1] + (-0.09960050292505122) * bodyOfConstraint.v_0[2] + 0.5648625214636235 * bodyOfConstraint.v_0[3]
3 : bodyOfConstraint.v_0[2] = 0.9663834860128886 * freeMotionScalarInit.v_rel_a_1 + 0.2370288965055821 * freeMotionScalarInit.v_rel_a_2 + (-0.09960050292505122) * freeMotionScalarInit.v_rel_a_3
4 : bodyOfConstraint.v_0[3] = (-0.1411940808771254) * freeMotionScalarInit.v_rel_a_1 + 0.8130157214783864 * freeMotionScalarInit.v_rel_a_2 + 0.5648625214636235 * freeMotionScalarInit.v_rel_a_3
[
  0.0 , -0.2148655148851041 , 0.0 , 1.0 ;
  -0.5648625214636235 , 0.0 , 0.09960050292505122 , -0.8191520442889918 ;
  0.0 , -0.9663834860128886 , 1.0 , 0.0 ;
  1.0 , 0.1411940808771254 , 0.0 , 0.0
]
  *
[
  bodyOfConstraint.v_0[3] ;
  freeMotionScalarInit.v_rel_a_1 ;
  bodyOfConstraint.v_0[2] ;
  bodyOfConstraint.v_0[1]
]
  =
[
  (-0.5318108111446446) * freeMotionScalarInit.v_rel_a_2 + 0.8191520442889918 * freeMotionScalarInit.v_rel_a_3 ;
  -0.0 ;
  0.2370288965055821 * freeMotionScalarInit.v_rel_a_2 + (-0.09960050292505122) * freeMotionScalarInit.v_rel_a_3 ;
  0.8130157214783864 * freeMotionScalarInit.v_rel_a_2 + 0.5648625214636235 * freeMotionScalarInit.v_rel_a_3
]
 might be structurally or numerically singular for variable bodyOfConstraint.v_0[1] since U(4,4) = 0.0. It might be hard to solve. Compilation continues anyway.

The LU-decomposition and determinant

L =

         0    0.2223    1.0000         0
   -0.5649   -0.0825   -0.8192    1.0000
         0    1.0000         0         0
    1.0000         0         0         0

U =

    1.0000    0.1412         0         0
         0   -0.9664    1.0000         0
         0         0   -0.2223    1.0000
         0         0         0   -0.0000
det =

   7.1565e-17

indeed show that this is close to singularity. I don't know what magic is applied but omc just solves this. Maybe the initial values are already solutions to this mess.

Note that ASSC only removes an equation when the LHS can be reduced to an actual zero. Even though we are working with Real (double) values and matlab computes a determinant that is almost but not quite zero, my algorithm lands the spot on zero. This is not close to but actual singularity.

When i apply ASSC on this, the algorithm correctly finds:

[ASSC] The equation: bodyOfConstraint.body.frame_a.r_0[3] = 0.2 + (-0.1411940808771254) * freeMotionScalarInit.derv[1].u + 0.8130157214783864 * freeMotionScalarInit.derv[2].u + 0.5648625214636235 * freeMotionScalarInit.derv[3].u
[ASSC] Gets replaced by equation: 0.0 = -0.3066830640280385 + 0.5648625214636235 * ((-0.1411940808771254) * (-0.2 + (-0.8191520442889918) * freeMotionScalarInit.derv[3].u + 0.5318108111446446 * freeMotionScalarInit.derv[2].u) + 0.2148655148851041 * (0.2 + 0.5648625214636235 * freeMotionScalarInit.derv[3].u - (-0.8130157214783864) * freeMotionScalarInit.derv[2].u)) + (-0.176326980708465) * (0.9663834860128886 * (-0.2 + (-0.8191520442889918) * freeMotionScalarInit.derv[3].u + 0.5318108111446446 * freeMotionScalarInit.derv[2].u) + 0.2148655148851041 * (-0.3 + (-0.09960050292505122) * freeMotionScalarInit.derv[3].u - (-0.2370288965055821) * freeMotionScalarInit.derv[2].u))

This creates a superfluous equation that gets removed by my new initialization. Instead an unfixed initialization variable is fixed:

removed eqns (1)
========================================
1/1 (1): 0.0 = -0.3066830640280385 + 0.5648625214636235 * ((-0.1411940808771254) * (-0.2 + (-0.8191520442889918) * freeMotionScalarInit.derv[3].u + 0.5318108111446446 * freeMotionScalarInit.derv[2].u) + 0.2148655148851041 * (0.2 + 0.5648625214636235 * freeMotionScalarInit.derv[3].u - (-0.8130157214783864) * freeMotionScalarInit.derv[2].u)) + (-0.176326980708465) * (0.9663834860128886 * (-0.2 + (-0.8191520442889918) * freeMotionScalarInit.derv[3].u + 0.5318108111446446 * freeMotionScalarInit.derv[2].u) + 0.2148655148851041 * (-0.3 + (-0.09960050292505122) * freeMotionScalarInit.derv[3].u - (-0.2370288965055821) * freeMotionScalarInit.derv[2].u))   [unknown |0|0|0|0|]

fixed vars (1)
========================================
1: bodyOfConstraint.body.frame_a.r_0[3]:VARIABLE(flow=false start = 0.0 unit = "m" fixed = true )  "Position vector from world frame to the connector frame origin, resolved in world frame" type: Real [3]

And as we can see the chosen variable makes sense. The singular algebraic loop has been ripped apart by removing one equation, so one of the previously involved variables has to be fixed.

Currently i have some follow up problems because ASSC seems to mess with some structures. As i said - probably ordering issues. So i cannot really test this right now but working on it. In theory i do not see any problem with my approach though.

@Francesco What is your opinion on this?

[Karim Adbdelhak]

I just realized that the singular system and my detected singularity are not actually the same.

The first one is a loop regarding the state derivatives v_0 and the second one ASSC detects is one regarding the states themselves r_0.

bodyOfConstraint.v_0[1] = der(bodyOfConstraint.body.frame_a.r_0[1])
bodyOfConstraint.v_0[2] = der(bodyOfConstraint.body.frame_a.r_0[2])
bodyOfConstraint.v_0[3] = der(bodyOfConstraint.body.frame_a.r_0[3])

The jacobians for both systems are the same though. So mathematically it is the same problem.

I am not sure whats happening but i will have a deeper look at this.

[Karim Adbdelhak]

Both are actually found:

[ASSC] The equation: 0.0 = 0.8191520442889918 * bodyOfConstraint.v_0[1] + (-0.09960050292505122) * bodyOfConstraint.v_0[2] + 0.5648625214636235 * bodyOfConstraint.v_0[3]                                                                              
[ASSC] Gets replaced by equation: 0.0 = 0.8191520442889918 * (0.9663834860128886 * ((-0.8191520442889918) * freeMotionScalarInit.v_rel_a_3 + 0.5318108111446446 * freeMotionScalarInit.v_rel_a_2) + 0.2148655148851041 * ((-0.09960050292505122) * free
MotionScalarInit.v_rel_a_3 - (-0.2370288965055821) * freeMotionScalarInit.v_rel_a_2)) + 0.5648625214636235 * ((-0.1411940808771254) * (0.09960050292505122 * freeMotionScalarInit.v_rel_a_3 + (-0.2370288965055821) * freeMotionScalarInit.v_rel_a_2) +
 0.9663834860128886 * (0.5648625214636235 * freeMotionScalarInit.v_rel_a_3 - (-0.8130157214783864) * freeMotionScalarInit.v_rel_a_2))

[Francesco Casella]

I had a look at the model, BTW, I attach a simplified version that replicates the problem with a lot less equations.

Dymola also finds two linear systems of four equations, that are solved symbolically via tearing:

  // Linear system of equations
  // Tag: initialization.linear[1]
    // Symbolic solution
      /*  Original equation
      (-5.551115123125783E-17)*freeMotionScalarInit.derv[1].u = (
        -0.036476928899629085)-(0.5318108111446445*freeMotionScalarInit.r_rel_a_2
        -0.8191520442889918*freeMotionScalarInit.r_rel_a_3+0.12158976299876358*(
        0.23702889650558212*freeMotionScalarInit.r_rel_a_2-0.0996005029250512*
        freeMotionScalarInit.r_rel_a_3-0.3)-0.6895698123465139*(0.8130157214783864
        *freeMotionScalarInit.r_rel_a_2+0.5648625214636234*freeMotionScalarInit.r_rel_a_3));
      */
      freeMotionScalarInit.derv[1].u := (-18014398509481984.0)*((
        -0.036476928899629085)-(0.5318108111446445*freeMotionScalarInit.r_rel_a_2
        -0.8191520442889918*freeMotionScalarInit.r_rel_a_3+0.12158976299876358*(
        0.23702889650558212*freeMotionScalarInit.r_rel_a_2-0.0996005029250512*
        freeMotionScalarInit.r_rel_a_3-0.3)-0.6895698123465139*(0.8130157214783864*
      freeMotionScalarInit.r_rel_a_2+0.5648625214636234*freeMotionScalarInit.r_rel_a_3)));

    // Torn part
      constraint.frame_b.r_0[3] := 0.2+0.14119408087712543*(7.082449871749603*(
        0.8130157214783864*freeMotionScalarInit.r_rel_a_2+0.5648625214636234*
        freeMotionScalarInit.r_rel_a_3)-freeMotionScalarInit.derv[1].u);
      constraint.frame_b.r_0[2] := 0.9663834860128886*freeMotionScalarInit.derv[1].u
        +0.23702889650558212*freeMotionScalarInit.r_rel_a_2-0.0996005029250512*
        freeMotionScalarInit.r_rel_a_3-0.3;
      constraint.frame_b.r_0[1] := (-1.220774588761456)*(0.5648625214636234*(
        constraint.frame_b.r_0[3]-0.2)-0.0996005029250512*constraint.frame_b.r_0[2]
        -0.19371055973531373);
  // End of linear system of equations


  // Linear system of equations
  // Tag: initialization.linear[2]
    // Symbolic solution
      /*  Original equation
      3.3306690738754696E-16*bodyOfConstraint.v_0[3] =  -(0.5318108111446445*
        freeMotionScalarInit.v_rel_a_2-0.8191520442889918*freeMotionScalarInit.v_rel_a_3
        +0.12158976299876358*(0.23702889650558212*freeMotionScalarInit.v_rel_a_2
        -0.0996005029250512*freeMotionScalarInit.v_rel_a_3)-0.6895698123465143*(
        0.8130157214783864*freeMotionScalarInit.v_rel_a_2+0.5648625214636234*
        freeMotionScalarInit.v_rel_a_3));
      */
      bodyOfConstraint.v_0[3] := (-3002399751580330.5)*(0.5318108111446445*
        freeMotionScalarInit.v_rel_a_2-0.8191520442889918*freeMotionScalarInit.v_rel_a_3
        +0.12158976299876358*(0.23702889650558212*freeMotionScalarInit.v_rel_a_2
        -0.0996005029250512*freeMotionScalarInit.v_rel_a_3)-0.6895698123465143*(
        0.8130157214783864*freeMotionScalarInit.v_rel_a_2+0.5648625214636234*
        freeMotionScalarInit.v_rel_a_3));

    // Torn part
      der(freeMotionScalarInit.derv[1].u) := 7.082449871749603*(0.8130157214783864
        *freeMotionScalarInit.v_rel_a_2+0.5648625214636234*freeMotionScalarInit.v_rel_a_3
        -bodyOfConstraint.v_0[3]);
      bodyOfConstraint.v_0[2] := 0.9663834860128886*der(freeMotionScalarInit.derv[1].u)
        +0.23702889650558212*freeMotionScalarInit.v_rel_a_2-0.0996005029250512*
        freeMotionScalarInit.v_rel_a_3;
      der(constraint.frame_b.r_0[1]) := (-1.220774588761456)*(0.5648625214636234
        *bodyOfConstraint.v_0[3]-0.0996005029250512*bodyOfConstraint.v_0[2]);
  // End of linear system of equations

As you can see, the coefficients multiplying the unknown of the residual equation are very, very small (-5.551115123125783E-17 and 3.3306690738754696E-16), but the system is solved anyway, and the solution makes sense. I guess that is what was also happening with omc.

I'm not an expert in this kind of models, so I can't really say why these initial equations turn out to be nearly singular. As I understand, after index reduction this model has four state variables, two positions and two velocities, related to the second and third component of the translational variable of the connectors, and the free-motion initializer has four initial conditions corresponding to them. So, I would expect the initialization system to be square, which is in fact the case, and well-posed, which unfortunately is not the case.

I opened ​#3875 on the MSL tracker, let's see if we can get some intel from there. However, the cases failing in Buildings are mostly thermo-fluid systems, which I understand best. @Karim, can you find out an example there?

[Francesco Casella]

Ah, one more thing, wouldn't it be safer if the ASSC algorithm can only combine equations with integer coefficients (e.g. the ones stemming from connection equations), so that we hava a guarantee that exact zero actually means what it says?

[Karim Adbdelhak]

Replying to casella:

Ah, one more thing, wouldn't it be safer if the ASSC algorithm can only combine equations with integer coefficients (e.g. the ones stemming from connection equations), so that we hava a guarantee that exact zero actually means what it says?

We can do that, i implemented it for real numbers at the end of 2020 because you requested it. :)

I will just hide it behind a flag so that the work is not lost. In some cases it might be beneficial.

[Francesco Casella]

Replying to Karim.Abdelhak:

We can do that, i implemented it for real numbers at the end of 2020 because you requested it. :)

Oops :)

Do you remember what was the rationale behind the request? I can't recall it now. In fact, it could be useful in case there are components that multiply or divide some flows by a number, multiplication would be an integer, but the division could be rendered as the multiplication by the inverse.

I will just hide it behind a flag so that the work is not lost. In some cases it might be beneficial.

Sounds good. I'm also eager to hear from the MultiBody experts about the case you brought up, I can't really make sense out of it.

[Karim Adbdelhak]

Replying to casella:

Do you remember what was the rationale behind the request? I can't recall it now.

The main rationale was that i still use actual zero comparison even for real valued coefficients (not with an epsilon, just flat on equality). The odd thing is that it triggered for this case. The only reason that comes to my mind is either another epsilon used in the generated C-code, or we hit our precision threshold and it is rounded in such a way that it actually ends up being zero.

Either way, i added the config flag --realASSC which defaults to false.

[Francesco Casella]

@Karim, it turned out that Modelica.Mechanics.MultiBody.Examples.Constraints.PrismaticConstraint actually had wrong singular initial conditions, see the discussion in MSL ​#7818. Apparently the problem is fixed on MSL master, can you check your algorithms with that version?

EDIT: These changes were also ported on the maint/3.2.3 and maint/4.0.1, see ​#7966, we will get them as soon as the package manager refers to them.

The question of course is what we should do in this cases. I guess the right solution is to warn the user of a nearly singular problem. If a solution is found, so be it, but one should be made aware that it is borderline. Dymola follows the approach of "least users's harassment" and remains mostly silent in these cases, but I don't think hiding potential problem is a good idea. At least from a researcher's perspective, I understand engineers in industry often need solutions that work, not silly questions the can't answer :)

The other question is: are the thermofluid examples failing in Buildings also plagued by bogus singular initial conditions? We should find out. If you can find some that has the same "nearly singular jacobian" problem, I can try to figure out what is the root cause.

[Francesco Casella]

Replying to Karim.Abdelhak:

Either way, i added the config flag --realASSC which defaults to false.

Was this harmful to the testuite?

[Karim Adbdelhak]

Replying to casella:

Replying to Karim.Abdelhak:

Either way, i added the config flag --realASSC which defaults to false.

Was this harmful to the testuite?

Not yet in the master, i will push it together with the other changes.

[Karim Adbdelhak]

Replying to Karim.Abdelhak:

Replying to casella:

Replying to Karim.Abdelhak:

Either way, i added the config flag --realASSC which defaults to false.

Was this harmful to the testuite?

Not yet in the master, i will push it together with the other changes. But it should not make any notable differences right now.

[Karim Adbdelhak]

Replying to casella:

The other question is: are the thermofluid examples failing in Buildings also plagued by bogus singular initial conditions? We should find out. If you can find some that has the same "nearly singular jacobian" problem, I can try to figure out what is the root cause.

For now i will try to push my new initialization, afterwards i will try for initial ASSC and following i will try to activate Real valued ASSC. Like this we can see what each part does step by step.

I will check with the current MSL version once my initialization is finalized, there is still some ordering issues that gives me headaches and i don't know why the compiler is so damn sensitive for that.

[Karim Adbdelhak]

I am struggling with some of the default examples in our testsuite for under/over/mixed determined initialization. Maybe you can enlighten me @Francesco

Take this one for example:

  model underdeterminedTest6
    Real a(start=-1);
    Real b(start=-2);

    Real x(start=-3);
    Real y(start=-4);

  initial equation 
    x = pre(a);
    y = pre(b);

  equation 
    der(x) = 0;
    der(y) = 0;

    when time > 0.5 then
      a = 2;
    end when;

    when initial() then
      b = 1;
    end when;
  end underdeterminedTest6;

The current initialization fixes x and y, but from my understanding that does not make any sense at all. Isn't it the case that discrete states (pre) are always set equal to their discrete variable counterparts during initialization (a=pre(a))? So in this case we would end up with an initial system of

  variables (8)
    a, b, x, y
    pre(a), pre(b), der(x), der(y)

  equations (7)
    a = pre(a);
    b = pre(b);
    x = pre(a);
    y = pre(b);
    der(x) = 0;
    der(y) = 0;
    when initial() then
      b = 1;
    end when;

A system of 8 variables and 7 equations. The only variables that could end up unmatched are either a or pre(a) depending on which way the first equation is solved. Since a is a suitable initialization variable and pre(a) is not, we would add the equation a = $START.a by fixing it.

It seems like the equations a = pre(a) are not currently forced upon the initialization system, but rather created if they have to be fixed. Is this correct and expected behaviour?

Does this model even make any sense?

[Francesco Casella]

Replying to Karim.Abdelhak:

I will check with the current MSL version once my initialization is finalized, there is still some ordering issues that gives me headaches and i don't know why the compiler is so damn sensitive for that.

Shouldn't we always order them lexicographically or something like that, so the compiler is insensitive to these issues?

[Francesco Casella]

Replying to Karim.Abdelhak:

The current initialization fixes x and y, but from my understanding that does not make any sense at all.

Of course it does :)

The system has four states, two continuous x and y and two discrete a and b. There are two initial equations, so you need two additional conditions, i.e. set fixed = true to two the four states.

Isn't it the case that discrete states (pre) are always set equal to their discrete variable counterparts during initialization (a=pre(a))?

Yes, but only if they're matched when equation is not active during initialization, see ​Section 8.6 of the spec.

The initial system is going to be

  variables (8)
    a, b, x, y
    pre(a), pre(b), der(x), der(y)
  equations (6)
    der(x) = 0;
    der(y) = 0;
    a = pre(a); // added because the when equation for a is not active @init
    b = 1;      // active @init because of if initial()
    x = pre(a);
    y = pre(b);

So you have two missing equations. The system can be made square by adding two initial equations involving x, y, pre(a), a, and pre(b). Fixing x and y looks good in order to make the initial system well-posed.

[Francesco Casella]

@Karim, regarding Modelica.Mechanics.MultiBody.Examples.Constraints.PrismaticConstraint, you wrote that you were using MSL 3.2, I'm not sure why you still consider that old version - it is not even Modelica compliant.

If you use the MSL 3.2.3 or 4.0.0 shipped with the package manager, you should get an updated version with the fixed non-singular initialization problem.

[Francesco Casella]

After ​PR #7903 was merged in, the three examples reported in this ticket

• ​Buildings.Fluid.HeatExchangers.Validation.HeaterCooler_u
• ​Buildings.Fluid.HeatExchangers.DXCoils.AirCooled.Examples.MultiStage
• ​Buildings.Experimental.DistrictHeatingCooling.Validation.HeatingCoolingHotWaterSmall

successfully pass the backend phase.