Unnecessary choice of dynamic state selection leads to code generation failure

[casella]

Please consider ​Modelica.Fluid.Examples.HeatingSystem.

The backend uses dynamic state selection

 * Number of states: 10 
($STATESET1.x[4],$STATESET1.x[3],$STATESET1.x[2],$STATESET1.x[1],
 radiator.mediums[1].p,radiator.mediums[1].h,pipe.mediums[1].p,
 pipe.mediums[1].h,pipe.mediums[2].p,pipe.mediums[2].h)

but then later fails to generate the code:

[/var/lib/jenkins/ws/OpenModelicaLibraryTestingWork/OpenModelica/OMCompiler/Compiler/SimCode/SimCodeUtil.mo:4450:9-4450:59:writable] Error: Internal error function createStateSetsSets failed.
[/var/lib/jenkins/ws/OpenModelicaLibraryTestingWork/OpenModelica/OMCompiler/Compiler/SimCode/SimCodeUtil.mo:737:5-737:146:writable] Error: Internal error function createSimCode failed [Transformation from optimised DAE to simulation code structure failed]

Dymola manages to make the following static state selection:

Statically selected continuous time states
heater.mediums[1].h
heater.mediums[1].p
pipe.mediums[1].h
pipe.mediums[1].p
pipe.mediums[2].h
pipe.mediums[2].p
radiator.mediums[1].h
radiator.mediums[1].p
tank.level
tank.medium.h

which is also the one I would expect as a modeller for this system. In fact, it reports that

We cannot differentiate equations for the following variables
having stateSelect = StateSelect.prefer.
They will be treated as if they had stateSelect = StateSelect.default:
pump.medium.h
pump.medium.p

which may be relevant to the discussion.

Karim, Andreas, Lennart, do you think you can fix the state selection algorithm so that it leads to the above-mentioned fixed state selection?

[casella]

Any progress on this issue? Should we reschedule it to 2.0.0?

[Karim.Abdelhak]

-d=bltdump

[Karim.Abdelhak]

Since the output of -d=bltdump is rather big, i added it as an additional file.

[casella]

@Karim, I will need you help to interpret the -d=bltdump output, we should organize a webmeeting next week, if possible, so you can explain me exactly how the state selection results from the reported structural analysis.

In the meantime, I noticed several odd things:

• Based on my understanding of the underlying physical modelling, I would expect that the only equations to be differentiated are related to mass and energy storage, so they may involve quantities such as mass, energy, density, specific enthalpy. That's what I see using Dymola. The set of equations in the OMC bltdump instead includes numerous equations describing mass and energy flows. I have no idea why this happens, but it doesn't feel right - mass and energy flows are not differentiated in this model, so I don't understand why they should end up in the state constraint equations.
• Some differentiated equations are clearly wrong, e.g.:

  ------------------53------------------
  Constraint equation to be differentiated:
  pump.eta = homotopy(Modelica.Fluid.Machines.ControlledPump$pump.efficiencyCharacteristic(pump.V_flow_single * pump.N_nominal / pump.N, 0.8), Modelica.Fluid.Machines.ControlledPump$pump.efficiencyCharacteristic(pump.V_flow_single_init, 0.8))
  Differentiated equation:
  der(pump.eta) = 0.0
  

• One differentiated equation even involves the port penetrations

  Constraint equation to be differentiated:
  tank.ports_penetration[2] = Modelica.Fluid.Utilities.regStep(tank.level + (-0.1) * tank.portsData[2].diameter - tank.portsData[2].height, 1.0, 0.001, 0.1 * tank.portsData[2].diameter)
  

  this also looks very odd to me.

It looks like we have some fundamental issues with the state selection algorithm - I think we should try to sort this thing out with high priority. Maybe this will solve other issues (e.g. #5328).

[casella]

Replying to casella:

• Some differentiated equations are clearly wrong, e.g.:

  ------------------53------------------
  Constraint equation to be differentiated:
  pump.eta = homotopy(Modelica.Fluid.Machines.ControlledPump$pump.efficiencyCharacteristic(pump.V_flow_single * pump.N_nominal / pump.N, 0.8), Modelica.Fluid.Machines.ControlledPump$pump.efficiencyCharacteristic(pump.V_flow_single_init, 0.8))
  

In fact, the function in question has a body that computes a constant efficiency output, so in this case OMC is actually smart to figure out that the derivative is indeed identically zero.

[casella]

The conclusion of today's discussion is that most likely the set of unmatched equations that backend outputs (42 of them) is too large, probably due to some issues that cause the loss of some dependencies in other equations, making the 42 currently selected equation orphans. In fact, the set of those equations is later proven to be unsolvable, which explains the failure of the state selection algorithm.

We should investigate what makes the matching go wrong, possibly on a simpler case that reproduces the issue. @casella will provide it ASAP.

[casella]

I attach a simplified version of HeatingSystem that only has four states. Dymola picks the following fixed state selection

Statically selected continuous time states
heater.mediums[1].h
heater.mediums[1].p
tank.level
tank.medium.h

while OMC selects

Number of states: 4 ($STATESET1.x[4],$STATESET1.x[3],$STATESET1.x[2],$STATESET1.x[1])

and later fails with

[12] 00:46:43 Translation Error
[C:/dev/OM64bit/OMCompiler/Compiler/SimCode/SimCodeUtil.mo: 4458:9-4458:59]: Internal error function createStateSetsSets failed.

[casella]

HeatingSystem2.mo is even smaller. Dymola selects:

Statically selected continuous time states
tank.level
tank.medium.h

while OMC outputs

Number of states: 2 ($STATESET1.x[2],$STATESET1.x[1])

then fails with

[13] 01:17:14 Translation Error
[C:/dev/OM64bit/OMCompiler/Compiler/SimCode/SimCodeUtil.mo: 4458:9-4458:59]: Internal error function createStateSetsSets failed.

Karim, I added -d=bltdump but didn't get any information. Can you please do that an attach the output to the ticket so I can have a look? Thanks!

[Karim.Abdelhak]

-d=bltdump of HeatingSystem2.mo

[Karim.Abdelhak]

Replying to casella:

Karim, I added -d=bltdump but didn't get any information.

Are you using the current nightly build? It should work.

Can you please do that an attach the output to the ticket so I can have a look? Thanks!

Yes. It looks like these are the equations that cause problems:

Unassigned equations: (3)
========================================
1/1 (1): 0.0 = if tank.regularFlow[1] then tank.ports[1].p - homotopy(tank.vessel_ps_static[1] + 0.5 * tank.portAreas[1] ^ (-2.0) * smooth(2, if noEvent((-pump.m_flow) >= tank.m_flow_small) then (-1.0 + tank.portsData[1].zeta_in + (tank.portAreas[1] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[1] / tank.portInDensities[1] * pump.m_flow ^ 2.0 else if noEvent((-pump.m_flow) <= (-tank.m_flow_small)) then pump.m_flow ^ 2.0 * (-1.0 + (tank.portAreas[1] / tank.vesselArea) ^ 2.0 - tank.portsData[1].zeta_out) / (tank.heatTransfer.states[1].d * tank.ports_penetration[1]) else if noEvent((-1.0 + tank.portsData[1].zeta_in + (tank.portAreas[1] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[1] / tank.portInDensities[1] >= (1.0 + tank.portsData[1].zeta_out - (tank.portAreas[1] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[1])) then Modelica.Fluid.Utilities.regSquare2.regSquare2_utility(-pump.m_flow, tank.m_flow_small, (-1.0 + tank.portsData[1].zeta_in + (tank.portAreas[1] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[1] / tank.portInDensities[1], (1.0 + tank.portsData[1].zeta_out - (tank.portAreas[1] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[1]), false, 1.0) else -Modelica.Fluid.Utilities.regSquare2.regSquare2_utility(pump.m_flow, tank.m_flow_small, (1.0 + tank.portsData[1].zeta_out - (tank.portAreas[1] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[1]), (-1.0 + tank.portsData[1].zeta_in + (tank.portAreas[1] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[1] / tank.portInDensities[1], false, 1.0)), tank.vessel_ps_static[1]) else if tank.inFlow[1] then tank.ports[1].p - tank.vessel_ps_static[1] else -pump.m_flow   [unknown |0|0|0|0|]
2/2 (1): 0.0 = if tank.regularFlow[2] then tank.ports[2].p - homotopy(tank.vessel_ps_static[2] + 0.5 * tank.portAreas[2] ^ (-2.0) * smooth(2, if noEvent(pump.m_flow >= tank.m_flow_small) then (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2] * pump.m_flow ^ 2.0 else if noEvent(pump.m_flow <= (-tank.m_flow_small)) then pump.m_flow ^ 2.0 * (-1.0 + (tank.portAreas[2] / tank.vesselArea) ^ 2.0 - tank.portsData[2].zeta_out) / (tank.heatTransfer.states[1].d * tank.ports_penetration[2]) else if noEvent((-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2] >= (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[2])) then Modelica.Fluid.Utilities.regSquare2.regSquare2_utility(pump.m_flow, tank.m_flow_small, (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2], (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[2]), false, 1.0) else -Modelica.Fluid.Utilities.regSquare2.regSquare2_utility(-pump.m_flow, tank.m_flow_small, (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[2]), (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2], false, 1.0)), tank.vessel_ps_static[2]) else if tank.inFlow[2] then tank.ports[2].p - tank.vessel_ps_static[2] else pump.m_flow   [unknown |0|0|0|0|]
3/3 (1): pump.Hb_flow = pump.m_flow * (smooth(0, tank.medium.h) - smooth(0, pump.medium.h))   [dynamic |0|0|0|0|]

They cannot be assigned statically and therefore dynamic state selection is necessary. As we discussed last thursday, they should not be considered constraint equations in the first place, right?

[casella]

Replying to Karim.Abdelhak:

Replying to casella:

Karim, I added -d=bltdump but didn't get any information.

Are you using the current nightly build?

No, sorry, I went back and forth with installed versions because of ticket:4504#comment:12, I guess I was using the "wrong" one. Each time uninstalling and reinstalling on Windows takes about 15 minutes, it's a bit of a pain.

Yes. It looks like these are the equations that cause problems:

Unassigned equations: (3)
========================================
1/1 (1): 0.0 = if tank.regularFlow[1] then tank.ports[1].p - homotopy(tank.vessel_ps_static[1] + 0.5 * tank.portAreas[1] ^ (-2.0) * smooth(2, if noEvent((-pump.m_flow) >= tank.m_flow_small) then (-1.0 + tank.portsData[1].zeta_in + (tank.portAreas[1] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[1] / tank.portInDensities[1] * pump.m_flow ^ 2.0 else if noEvent((-pump.m_flow) <= (-tank.m_flow_small)) then pump.m_flow ^ 2.0 * (-1.0 + (tank.portAreas[1] / tank.vesselArea) ^ 2.0 - tank.portsData[1].zeta_out) / (tank.heatTransfer.states[1].d * tank.ports_penetration[1]) else if noEvent((-1.0 + tank.portsData[1].zeta_in + (tank.portAreas[1] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[1] / tank.portInDensities[1] >= (1.0 + tank.portsData[1].zeta_out - (tank.portAreas[1] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[1])) then Modelica.Fluid.Utilities.regSquare2.regSquare2_utility(-pump.m_flow, tank.m_flow_small, (-1.0 + tank.portsData[1].zeta_in + (tank.portAreas[1] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[1] / tank.portInDensities[1], (1.0 + tank.portsData[1].zeta_out - (tank.portAreas[1] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[1]), false, 1.0) else -Modelica.Fluid.Utilities.regSquare2.regSquare2_utility(pump.m_flow, tank.m_flow_small, (1.0 + tank.portsData[1].zeta_out - (tank.portAreas[1] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[1]), (-1.0 + tank.portsData[1].zeta_in + (tank.portAreas[1] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[1] / tank.portInDensities[1], false, 1.0)), tank.vessel_ps_static[1]) else if tank.inFlow[1] then tank.ports[1].p - tank.vessel_ps_static[1] else -pump.m_flow   [unknown |0|0|0|0|]
2/2 (1): 0.0 = if tank.regularFlow[2] then tank.ports[2].p - homotopy(tank.vessel_ps_static[2] + 0.5 * tank.portAreas[2] ^ (-2.0) * smooth(2, if noEvent(pump.m_flow >= tank.m_flow_small) then (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2] * pump.m_flow ^ 2.0 else if noEvent(pump.m_flow <= (-tank.m_flow_small)) then pump.m_flow ^ 2.0 * (-1.0 + (tank.portAreas[2] / tank.vesselArea) ^ 2.0 - tank.portsData[2].zeta_out) / (tank.heatTransfer.states[1].d * tank.ports_penetration[2]) else if noEvent((-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2] >= (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[2])) then Modelica.Fluid.Utilities.regSquare2.regSquare2_utility(pump.m_flow, tank.m_flow_small, (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2], (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[2]), false, 1.0) else -Modelica.Fluid.Utilities.regSquare2.regSquare2_utility(-pump.m_flow, tank.m_flow_small, (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[2]), (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2], false, 1.0)), tank.vessel_ps_static[2]) else if tank.inFlow[2] then tank.ports[2].p - tank.vessel_ps_static[2] else pump.m_flow   [unknown |0|0|0|0|]
3/3 (1): pump.Hb_flow = pump.m_flow * (smooth(0, tank.medium.h) - smooth(0, pump.medium.h))   [dynamic |0|0|0|0|]

They cannot be assigned statically and therefore dynamic state selection is necessary. As we discussed last thursday, they should not be considered constraint equations in the first place, right?

Absolutely.

[casella]

I think I finally understood what the problem is.

I picked the relevant equations of the system, which are the mass and energy balance equations of the pump and of the valve. Notice that the former are dynamic, while the latter are static, because of the choice of structural parameters:

der(tank.m) = tank.mb_flow;
der(tank.U) = tank.Hb_flow + tank.Qb_flow + tank.Wb_flow;
tank.m = tank.fluidVolume * tank.medium.d;U = m*medium.u;
tank.medium.u = tank.medium.h - tank.medium.p / tank.medium.d
tank.fluidVolume = tank.V; (binding equation)
tank.V = tank.crossArea * tank.level;
tank.medium.p = p_ambient;

0.0 = pump.Hb_flow + pump.Qb_flow + pump.Wb_flow;
0.0 = pump.mb_flow;
pump.fluidVolume = 0.0 (binding equation)
pump.m = pump.fluidVolume * pump.medium.d;
pump.U = pump.m * pump.medium.u;
pump.medium.u = pump.medium.h - pump.medium.p / pump.medium.d;

Now, stateSelect is set to StateSelect.prefer for the following variables:

tank.level
tank.medium.p
tank.medium.h
pump.medium.p
pump.medium.h

In fact, the pump is set to not have any internal mass and energy storage, hence there are no der(pump.m) nor der(pump.U) terms in the mass and energy equations. In fact pump.m = 0 and pump.U = 0.

In this case it makes obviously no sense to try to use pump.medium.p and pump.medium.h as states; pump.medium.preferredMediumStates should have been set to false, so that the last two variables do not have StateSelect.prefer, but rather StateSelect.default. Anyway, if a static state selection is possible, I guess the tool should pick it even if this implies not selecting some states which have stateSelect.prefer.

Regarding the tank, stateSelect.prefer was added to tank.level, without bothering of setting stateSelect.default on tank.medium.p, which is made trivial by the equation tank.medium.p = p_ambient;

It looks like the back-end is instead trying to push a bit too hard to select all of these variables with stateSelect.prefer as states.

In fact, if I add the modifier medium(preferredMediumStates = false) to pump (see the attached HeatingSystem3), the backend makes the correct state selection

Number of states: 2 (tank.level,tank.medium.h)

and the system is simulated correctly.

Again, one may argue that this is a modelling problem. However, the StateSelect.prefer attribute is only a hint to prefer the variable, but it is not an obligation to use it. It seems to me that the current state selection algorithm gives too much weigth to this setting, in those cases in which there is some redundancy.

@Karim, can you elaborate a bit further on how the stateSelect attribute is used in the formation of the MSSS and in the state selection algorithm?

I hope this helps you to fix the problem. As you and @lochel already guessed, I am also convinced that this issue is present in a number of other models in the MSL and in other libraries, so fixing it should have a beneficial effect. Apparently, the state selection algorithm in Dymola is more robust, so people inadvertently set redundant stateSelect.prefer attributes without realizing it, because Dymola does not complain and manages the models correctly. I guess this happens in many models out there, so we should try to support this scenario.

If necessary, I'm available this week for a further round of discussion.

[casella]

Andrea has a nice test case in which the state selection algorithm is using the information provided by the StateSelect attributes in a very questionable way. We're currently trying to reduce it to the barest minimum number of equations (could be as few as 3 or 4 DAEs).

I think it could be very useful to debug the state selection algorithm, we'll post it later today. Stay tuned.

[casella]

Please check the attached TestStateSelect package. The test model has one differential equation, a trivial dynamic energy balance

N*MM*Cp*der(T) - Vg*der(p)=0;

with two differentiated variables, p and T, and one constraint equation between them, which is the Van der Waals equation assuming constant number of moles N

N*R*T = (p + a*N^2/Vg^2) * (Vg - N*b);

The models in sub-package UnknownParameters compute the actual value of N via an initial equation, while the models in the other sub-package get the value from a binding equation.

Looking in the UnknownParameters package, we have four cases.

1. Pdefault_Tdefault has StateSelect.default on both differentiated variables. OMC chooses a dynamic state selection set with both p and T. It is true that variable Vg could become zero, so a fixed state selection on p may not be always a good idea for the compiler. In fact, Vg will never get to be zero, but the compiler doesn't know it, unless it goes through some non-trivial symbolic analysis. On the other hand, Van der Waals' equation can be solved for p assuming Vg > N*b, but the compiler doesn't know about this. So, a dynamic state selection seems reasonable.

2. Pprefer_Tavoid: in this case we set StateSelect.prefer on p and StateSelect.avoid on T. As I understand it from a modeller's point of view, this means: I know what I'm doing (-> Vg will never be zero), please use p as a state and make T a dummy state. Unfortunately, OMC still goes for a dynamic state selection. From my point of view this is questionable, but on the other hand that's also what Dymola does (see below).

3. Pprefer_Tnever: here we try to tell OMC really not to use T as a state variable, but we still get it in the dynamic state selection. This is definitely wrong, never means never.

4. Palways_Tnever: only by setting StateSelect.always on p we are actually able to convince OMC to make T a dummy state and go for a fixed state selection on p. However, we get the most puzzling warning message: " Variable $DER.p has attribute stateSelect=StateSelect.always, but was selected as a continuous variable.". Why but?? This really makes no sense at all and smells of bug from very far.

If we now instead give the value of N in a binding equation (package KnownParameter), in all four cases we get a static state selection on p. This is in fact good from a practical point of view, but I don't understand what is the difference between having the value of the parameter computed at runtime by initial equations, or having it provided by the user in the xml init file. From my point of view, in terms of structural analysis of the dynamic equation, there is no difference at all. So, I find it weird that the state selection algorithm produces different results based on this fact.

Summing up, I would expect to see the following results

• Pdefault_Tdefault: selecting p is computationally the best choice, except when Vg becomes zero; selecting T requires to solve Van der Waals' equation for p, which is only possible assuming a priori that Vg - N*b > 0, which the compiler can't be sure about, unless there is an assert somewhere. I would say dynamic state selection is the best option here, since it can use p most of the time, which is easier to solve but still be safe against Vg becoming zero.
• Pprefer_Tavoid: here we should clearly get p statically selected
• Pprefer_Tnever: here we should clearly get p statically selected
• Palways_Tnever: here we should clearly get p statically selected, without weird warning messages

The outcome from the two sub-packages should be the same, IMHO there is no reason why it should be different, as N could become zero with equal likelyhood in both cases.

I also tried these test cases with Dymola 2020. These are the results

• UnknownParameter
  • Pdefault_Tdefault: dynamic selection (ok)
  • Pprefer_Tavoid: dynamic selection (this is questionable, but maybe it is a matter of interpretation how strong prefer and avoid are)
  • Pprefer_Tnever: static selection p (ok)
  • Palways_Tnever: static selection p (ok
• KnownParameter
  • Pdefault_Tdefault: static selection p
  • Pprefer_Tavoid: static selection p
  • Pprefer_Tnever: static selection p
  • Palways_Tnever: static selection p

Also Dymola treats the two cases differently - maybe there is some good reason to do so, but honestly I can't see that.

[Karim.Abdelhak]

Replying to casella:

It seems like the stateSelection attribute is not used correctly, i will look into that. But i am not quite sure about the following point you made:

• Pprefer_Tavoid: here we should clearly get p statically selected

If we handle it like this there would be no difference between prefer and always and no difference between avoid and never. From what i understand prefer and avoid only influence the behaviour in the case of static selection. I there is a choice - prefer the one with the better attribute. Static selection can only be (locally) enforced by always and general selection prevented by never.

Are you sure Pprefer_Tavoid should lead to static selection? I think dymola is right here.

[Karim.Abdelhak]

OMC has a heuristic behind it which rates each state. The stateSelection attribute has following influence:

          case DAE.NEVER() then -20.0;
          case DAE.AVOID() then -1.5;
          case DAE.DEFAULT() then 0.0;
          case DAE.PREFER() then 1.5;
          case DAE.ALWAYS() then 20.0;

There are some special cases for DAE.ALWAYS(), where states get statically chosen, but none for DAE.NEVER() such that dynamic state selection can be prevented.

Some information on the selection heuristic can be looked at with -d=dummyselect. The result for Pprefer_Tnever is:

Prio 1 : 0.0
Prio 2 : 0.0
Prio 3 : 0.0
Prio 4 : 0.0
Prio 5 : 0.075
Calc Prio for Vg
 Prio StateSelect : 0.0
 Prio Heuristik : 0.075
 ### Prio Result : 0.075
Prio 1 : 0.1
Prio 2 : 0.0
Prio 3 : 0.0
Prio 4 : 0.0
Prio 5 : 0.0
Calc Prio for p
 Prio StateSelect : 1.5
 Prio Heuristik : 0.1
 ### Prio Result : 1.6
Prio 1 : 0.1
Prio 2 : 0.0
Prio 3 : 0.0
Prio 4 : 0.0
Prio 5 : 0.0
Calc Prio for T
 Prio StateSelect : -20.0
 Prio Heuristik : 0.1
 ### Prio Result : -19.9

[Karim.Abdelhak]

I found something you might not be aware of which can already be seen in the -d=dummyselect dump i presented. Vg also is considered a state candidate.

I will attach a file with the informations of -d=bltdump where you can see that following equation gets derived:

------------------1------------------
Constraint equation to be differentiated:
Vg = 0.125 + 0.01 * sin(6.283 * time)
Differentiated equation:
der(Vg) = 0.06283000000000001 * cos(6.283 * time)

This results in Vg being a state. Preventing T from being chosen still does not suffice since omc can still chose between p and Vg. Somehow dymola seems to keep those to stateSets apart and we don't.
Further information in attached file.

[Karim.Abdelhak]

I actually made an update that fixes this particular problem, but broke two other examples. I looked deeper into it and saw that N is handled as a regular variable rather than a parameter in the backend.

I will investigate if it is a frontend, backend or a transition problem. In the model it is clearly stated as final parameter but i guess since it can't be evaluated before runtime it is handled differently.

[casella]

Replying to Karim.Abdelhak:

It seems like the stateSelection attribute is not used correctly, i will look into that. But i am not quite sure about the following point you made:

• Pprefer_Tavoid: here we should clearly get p statically selected

If we handle it like this there would be no difference between prefer and always and no difference between avoid and never.

This is the definition of the StateSelect attributes in the language specification

  type StateSelect = enumeration(
    never "Do not use as state at all.",
    avoid "Use as state, if it cannot be avoided (but only if variable appears
           differentiated and no other potential state with attribute   
           default, prefer, or always can be selected).",
    default "Use as state if appropriate, but only if variable appears
             differentiated.",
    prefer "Prefer it as state over those having the default value
            (also variables can be selected, which do not appear
             differentiated). ",
    always "Do use it as a state."
);

Note that the specification does not explicitly mention the concept of static vs. dynamic state selection - this is left to the tool implementors.

The way I always interpreted this (but I was possibly wrong, given how Dymola acts), is that always means you must select this variable as state, so that if you can't do that, you end up with an error. So

model M
  Real x(StateSelect.always);
  Real y(StateSelect.always);
  Real q;
equation
  der(x) = -q + 1;
  der(y) = q;
  x = y;
end M;

is invalid, because you can't have both x and y as states at the same time. From this point of view, I think that having a dynamic state selection with one state chosen between x and y would just be cheating. Conversely, if you used prefer, then it would be ok to select one among x and y statically as state variable.

For this reason, I am always hesitant at selecting StateSelect.always on a pressure variable in a thermo-hydraulic model, because if this pressure ends up being equal to the pressure on the connector, you cannot be 100% sure that some other component will be connected to the connector that also wants to always have the (same) pressure as state. So, I basically only use always for debugging: if I want a certain variable to be a state, but I can't, then I force it with always, and if I get a failure I can check the reason for the failure and get some diagnostic information - for example, there could be a wrong equation somewhere that actually prevents p to be a state, and I can find that out.

So, my interpretation of prefer and always is that with prefer you try hard to select that as a state, possibly taking some chances (e.g. a possible division by zero in some cases), but if you can't do that for structural reasons (see the example above), then it's ok not to have it as a state; with always, in that case the compilation fails.

Are you sure Pprefer_Tavoid should lead to static selection? I think dymola is right here.

Following the specification to the letter, p should be preferred as a state over those having the default value (I guess this holds a fortiori if they have avoid), while T should be used as state, if it cannot be avoided, but only if no other potential state with attribute default, prefer, or always can be selected. I understand this clearly rules out a dynamic state selection where p and T are on the same level, the only way to enforce these requirements is to statically select p. Do I miss something?

Maybe we should discuss this on the modelica-design list. I'm sure it will be a long thread :)

[casella]

Replying to Karim.Abdelhak:

OMC has a heuristic behind it which rates each state. The stateSelection attribute has following influence:

          case DAE.NEVER() then -20.0;
          case DAE.AVOID() then -1.5;
          case DAE.DEFAULT() then 0.0;
          case DAE.PREFER() then 1.5;
          case DAE.ALWAYS() then 20.0;

I'm not sure how these weights are actually used, but I guess NEVER should be much lower and PREFER much higher, otherwise there is a finite chance of them getting selected (or not selected) against the attribute.

Maybe PREFER and AVOID should have higher weights, so they are not selected (or selected) only if there are very strong reasons not to do so?

There are some special cases for DAE.ALWAYS(), where states get statically chosen

Good.

but none for DAE.NEVER() such that dynamic state selection can be prevented.

That's not what the specification says: "Do not use as state at all". This clearly prevents any exception. In German you'd say "Verboten" and that would be it :)

[casella]

Replying to Karim.Abdelhak:

I found something you might not be aware of which can already be seen in the -d=dummyselect dump i presented. Vg also is considered a state candidate.

This is quite odd. I mean, Vg is completely determined as a function of time.

Vg = 0.125 + 0.01 * sin(6.283 * time)

How on earth could it be a state?

Besides, I don't get why this equation should be unmatched, and thus end up in the MSSS sets. We have one equation with just one unknown Vg. The matching is very obvious to me. Do I miss something?

I will attach a file with the informations of -d=bltdump where you can see that following equation gets derived:

------------------1------------------
Constraint equation to be differentiated:
Vg = 0.125 + 0.01 * sin(6.283 * time)
Differentiated equation:
der(Vg) = 0.06283000000000001 * cos(6.283 * time)

This results in Vg being a state. Preventing T from being chosen still does not suffice since omc can still chose between p and Vg. Somehow dymola seems to keep those to stateSets apart and we don't.

I still don't get it. Pantelides' algorithm was invented to answer the following question: "what is the maximum number of initial values that I can arbitrarily choose in a DAE, and on which variables?". It is clear that Vg cannot be assigned any arbitrary initial value, so why should it be ever considered as a state candidate?

[casella]

Replying to Karim.Abdelhak:

I actually made an update that fixes this particular problem, but broke two other examples. I looked deeper into it and saw that N is handled as a regular variable rather than a parameter in the backend.

That explains the observed behaviour.

I will investigate if it is a frontend, backend or a transition problem. In the model it is clearly stated as final parameter but i guess since it can't be evaluated before runtime it is handled differently.

In fact, all non-structural parameters share the same property, so I don't see why it should.

[casella]

See also #4483, may be somehow related

[Karim.Abdelhak]

Replying to casella:

This is quite odd. I mean, Vg is completely determined as a function of time.

Vg = 0.125 + 0.01 * sin(6.283 * time)

How on earth could it be a state?

When looking at the bipartite graph

you can actually see why omc derives the equation. Since you have to derive equation (2) there will be a derivative of Vg in the system, resulting in it being a potential state. Right afterwards it gets chosen as dummy derivative due to the equation you mentioned, therefore it is just the case that an equation for der(Vg), which is actually also just an algebraic variable, is needed.

The MSSS here is (1,2) - (Vg). This seems correct to me.

[Karim.Abdelhak]

I actually just found another problem we have with these models.

The equation

 N*R*T = (p + a*N^2/Vg^2) * (Vg - N*b);

does not seem to get checked for the RHS. T could not be chosen because N was not detected to be a parameter. I fixed that (will be pushed soon i hope) in combination with states having stateSelect.never being forced as dummy derivatives. Now even in the default/default model states get statically selected since the RHS does not get checked if it devides by a variable in the first place. T should not be selected, since the RHS devides by Vg.

I don't quite know if this is only the case for state selection or a general problem in matching algorithms we have, but maybe this is worth a separate ticket. Once again... going to look deeper into it.

EDIT: p also can't be solved in this equation without deriving by Vg. It's even worse since (Vg - N*b) also is a divisor right?. I am not quite sure how this should be handled, maybe we need a new concept for these things.

[casella]

Replying to Karim.Abdelhak:

When looking at the bipartite graph

you can actually see why omc derives the equation. Since you have to derive equation (2) there will be a derivative of Vg in the system, resulting in it being a potential state. Right afterwards it gets chosen as dummy derivative due to the equation you mentioned, therefore it is just the case that an equation for der(Vg), which is actually also just an algebraic variable, is needed.

You're right, as long as der(Vg) pops up in the diffentiated equations, it has to be a potential state. But then, it is immediately removed as dummy.

The MSSS here is (1,2) - (Vg). This seems correct to me.

You're right again, I misinterpreted the log output.

[casella]

Replying to Karim.Abdelhak:

I actually just found another problem we have with these models.

The equation

 N*R*T = (p + a*N^2/Vg^2) * (Vg - N*b);

does not seem to get checked for the RHS.

Checked against what, exactly?

T could not be chosen because N was not detected to be a parameter. I fixed that (will be pushed soon i hope) in combination with states having stateSelect.never being forced as dummy derivatives.

Sounds good.

Now even in the default/default model states get statically selected

I understand p is selected, right?

since the RHS does not get checked if it devides by a variable in the first place. T should not be selected, since the RHS devides by Vg.

I don't quite know if this is only the case for state selection or a general problem in matching algorithms we have, but maybe this is worth a separate ticket. Once again... going to look deeper into it.

EDIT: p also can't be solved in this equation without deriving by Vg. It's even worse since (Vg - N*b) also is a divisor right?. I am not quite sure how this should be handled, maybe we need a new concept for these things.

I'm not sure I'm following you completely. If p is statically selected as a state, then you need to solve the equation for T, resulting in

T := 1/(N*R)*(p + a*N^2/Vg^2)*(Vg - N*b);

This may cause division by zero if the parameter N is later set to zero, but that's a chance we always take in OMC. For example

C*der(T) = Q;

will always be solved as

der(T) = 1/C*Q

and if C is set to zero, you get a division by zero, and that's too bad for you. If you want to handle that case, you need to evaluate the parameter, so that when C = 0, der(T) is symbolically eliminated from the system.

The fact that the RHS of the original equation (not of the solved one) has a division by Vg^2 that could potentially throw a division by zero exception at runtime is a fact of life. Van der Waal's equation is valid in a certain range of values of Vg, and if you go outside that range, the model will break. No amount of dynamic state selection or other wizardry can change that. Note that in this case the division was written by the modeller, so he/she is aware of the potential issues and is fine with them. In fact, he/she could add some min/max attributes to prevent the division by zero and generate a meaningful exception before that happens, but that's another story.

If you select T as state, instead, you need to solve the equation for p. This can be done in closed form:

p = N*R*T/(Vg - N*b)- (p + a*N^2/Vg^2)

this requires you to divide by the time-dependent expression (Vg - N*b), which is a hazard the compiler should probably not take. Note that this division is not present in the original equation.

If you couldn't solve this equation symbolically, you may miss this fact, unless you can determine that the Jacobian of the residual w.r.t. T} could become singular - normally this kind of check is not performed. In this case, the modeller may help preventing this potential problem by setting the StateSelect attribute to steer the backend towards selecting p, thus avoiding potential problems.

Anyway, in this case I would expect p to be selected if N is a parameter, in all cases, particularly as soon as p has prefer and/or T has avoid.

[casella]

Some further analysis during today's webmeeting: the unsolvable equation number 36 in the first log is expected to be implicitly solved after static state selection is made (there is no way to avoid that), so in the matching analysis we should also allow for variables that must be implicitly solved. Maybe that would solve the issue, it could possibly break other models.

[casella]

One more interesting piece of information would be to check which models in the ​library testsuite actually end up using dynamic state selection. From my understanding, only multibody systems with no direct connection to the fixed world frame should, I am not aware of any other libraries that actually need that. Maybe there indeed are some, but we should find out whether this is really necessary, on a case-by-case basis.

@adrpo, can you please add the -d=stateselection flag to the testsuite, and run grep -e="$STATESET" at the end of the test run, with output on the console. I can then fetch the console output from Jenkins and post-process it myself.

[casella]

See also #5574

[casella]

@Karim, before you leave for your holidays, could you summarize a plan to resolve this issue (if you have any?). I'm afraid that in a few weeks this will become too cold and we'll have to restart more or less from scratch.

Thanks!

[Karim.Abdelhak]

Replying to casella:

@Karim, before you leave for your holidays, could you summarize a plan to resolve this issue (if you have any?). I'm afraid that in a few weeks this will become too cold and we'll have to restart more or less from scratch.

Thanks!

I will only be on vacation for two weeks from tomorrow.

The basic idea was too look if all the stateSelect attributes get parsed correctly and based on that correct choices are made. I already pushed an update regarding StateSelect.never, but there may be more to do. Furthermore i will look at the heuristic which sorts the variables such that higher priority ones are more likely to be chosen and question if this is the way we want to do it. Maybe finding assignments of variables with different StateSelect attributes should be completely seperated such that it is ensured that the priorities are used correctly (We have ContinueMatching for that, why not use it?).

If this does not resolve the problem i might have to question the generation of the adjacency matrix in general. Just as in #5170 it might be the case that we don't inline function calls which have dummy inputs that are not actually used in the function at all. They should not be connected to the function in the bipartite graph because they can't be solved by it.

@Francesco, could you check the functions of the model if there are any unused inputs or stuff that should not end up in the adjacency matrix? Maybe there are also other cases i do not consider (record constructors etc.?).

[casella]

Replying to Karim.Abdelhak:

@Francesco, could you check the functions of the model if there are any unused inputs or stuff that should not end up in the adjacency matrix? Maybe there are also other cases i do not consider (record constructors etc.?).

Whoa, there are 122 function calls and constructors in the simples HeatingSystem3 model, going through them one by one is a bit inconvenient. Besides, they involve the Modelica.Fluid IF97 water model, which has a very tricky hidden caching systems based on LateInline and CSE. I wonder if that could also be part of the problem.

I searched for the keyword dummy in the flat model and I found that there are several functions with inputs that may or may not be used, such as viscosity, but my understanding is that no sane model should try to solve for the viscosity given other quantities, normally you compute the viscosity given pressure and temperature/specific enthalpy.

I think we should go once more through the bltdump and try to figure out why certain dummy candidate states are selected.

[casella]

BTW, I just retried with the latest nightly build. The simplified HeatingSystem2 and HeatingSystem3, which previously failed because of bogus dynamic state selection, now lead to static state selection, and successful simulation.

The simplified HeatingSystem attached to this ticket now has a weird behaviour. It first states up front

[2] 23:53:53 Symbolic Error
[Modelica.Fluid.Machines: 458:5-459:59]:
Model is structurally singular, error found sorting equations

then it continues (why? shouldn't the compilation abort after finding an error?), then it states again

[7] 23:53:53 Symbolic Error
Model is structurally singular, error found sorting equations 

but on a different set of equations, then it states:

[8] 23:53:53 Translation Error
Internal error Transformation Module PFPlusExt index Reduction Method Pantelides failed!

and reports a completely dynamic state selection

* Number of states: 4 ($STATESET1.x[4],$STATESET1.x[3],$STATESET1.x[2],$STATESET1.x[1])

finally failing with

[12] 23:53:53 Translation Error
[C:/dev/OM64bit/OMCompiler/Compiler/SimCode/SimCodeUtil.mo: 4459:9-4459:59]: Internal error function createStateSetsSets failed.

[13] 23:53:53 Translation Error
[C:/dev/OM64bit/OMCompiler/Compiler/SimCode/SimCodeUtil.mo: 746:5-746:146]: Internal error function createSimCode failed [Transformation from optimised DAE to simulation code structure failed]

I think we should investigate the reason why the backend initially reports the model to be singular. Up to my understanding, it is not; Dymola simulates it without any problems, using a fixed state selection. Maybe understanding this first bogus reported singularity holds the key to solving the riddle.

[Karim.Abdelhak]

I found the error. It is so shockingly simple that i did not consider it at first.

When parsing the frontend structure to backend structure, variables with StateSelect.prefer got directly created as BackendDAE.STATE(), which is clearly wrong. It should never be possible to create a BackendDAE.STATE() without knowing that the derivative actually exists in the system. That is why we ended up with all those weird equations in state selection.

An exception to this is, of course, StateSelect.always. It directly creates the variable as a BackendDAE.STATE() and the compilation fails if stateselection detects that it cannot be forced in the way the user intended to. This is expected behaviour i guess.

I updated it in ​PR381, the original HeatingSystem works with this commit (although homothopy is needed), let's see what jenkins has to say about this.

[Karim.Abdelhak]

I got 18 errors on the testsuite (​jenkins run), of which the following seem to be relevant:

1. The following models only break in MSL 3.1 and work in the current one. It seems like some min - max constraints don't work out the way they should. When running MixtureGases i realized that medium2.T equals medium2.T_degC which seems wrong, since T should be Kelvin and T_degC should be degrees Celsius. Do you know what exactly changed here from 3.1 to now?

   Modelica.Media.Examples.MixtureGases.mos – simulation_libraries_msl31_media
   Modelica.Media.Examples.TestOnly.IdealGasN2Mix.mos – simulation_libraries_msl31_media
   Modelica.Media.Examples.TestOnly.IdealGasN2.mos – simulation_libraries_msl31_media
   

2. This one breaks in the current one, it seems like the selection of a StateSelect.never state breaks the initialization.

   Modelica.Fluid.Examples.PumpingSystem.mos – simulation_libraries_msl32
   

3. Additional output from the NLS-solver, i am not quite sure about this one, since it returns the correct values.

   Modelica.Media.Examples.Tests.MediaTestModels.Water.WaterIF97OnePhase_ph.mos – simulation_libraries_msl32
   

4. Warnings about selecting StateSelect.never variables.

   Modelica.Fluid.Examples.Tanks.TanksWithOverflow.mos – simulation_libraries_msl32
   Modelica.Fluid.Examples.Explanatory.MeasuringTemperature.mos – simulation_libraries_msl32
   Modelica.Fluid.Examples.AST_BatchPlant.BatchPlant_StandardWater.mos – simulation_libraries_msl32
   

5. Optimization produces different linearization

   testSteamPipe.mos – openmodelica_linearization
   testDrumBoiler.mos – openmodelica_linearization
   

6. Similar to 1. an assert due to min/max gets triggered

   Ticket4258a.mos – simulation_modelica_inheritances
   

I guess the others are false negatives.

@Francesco Could you give me some thoughts on what i reported, especially to 1.? And is it possible that the selection of StateSelect.never states is actually correct or is it always wrong since there was a selection we could make without choosing them?

[Karim.Abdelhak]

I need help on one big problem. I updated a lot and most of it works (or i know how to make it work) but one mayor issue remains. I updated the handling of StateSelect.never vars so that they are forced in, but i got issues when running Modelica.Mechanics.MultiBody.Examples.Constraints.RevoluteConstraint. The following variables only appear linear dependent and the factors are time dependent so we cannot be sure that they become zero at some point:

StateSelect.never variables that will be forced as dummys (2)
========================================
1: freeMotionScalarInit.angle_d_1:STATE(1)(start = 0.0 unit = "rad/s" stateSelect=StateSelect.never )  "= der(     angle_1)" type: Real
2: freeMotionScalarInit.angle_d_2:STATE(1)(start = 0.0 unit = "rad/s" fixed = true stateSelect=StateSelect.nev     er )  "= der(angle_2)" type: Real

When forcing them in, all linear solvers fail at time 8.759925 (there are more confusing error messages but this is the most relevant one). I tried some different possibilities regarding pre-ordering equations and variables before matching to solve these, but none worked out.

Equations where freeMotionScalarInit.angle_d_1 appears:

16/16 (1): $DER.freeMotionScalarInit.initAngle.R_rel.T[2,3] = cos(freeMotionScalarInit.derd[3].u) * cos(freeMotionScalarInit.derd[1].u) * freeMotionScalarInit.angle_d_1 + sin(freeMotionScalarInit.derd[3].u) * (freeMotionScalarInit.initAngle.R_rel.T[3,1] * (-sin(freeMotionScalarInit.derd[1].u)) * freeMotionScalarInit.angle_d_1 + $DER.freeMotionScalarInit.initAngle.R_rel.T[3,1] * cos(freeMotionScalarInit.derd[1].u)) + cos(freeMotionScalarInit.derd[3].u) * freeMotionScalarInit.angle_d_3 * freeMotionScalarInit.initAngle.R_rel.T[3,1] * cos(freeMotionScalarInit.derd[1].u) - sin(freeMotionScalarInit.derd[3].u) * freeMotionScalarInit.angle_d_3 * sin(freeMotionScalarInit.derd[1].u)   [dynamic |0|0|0|0|]
17/17 (1): $DER.freeMotionScalarInit.initAngle.R_rel.T[2,2] = cos(freeMotionScalarInit.derd[3].u) * (-sin(freeMotionScalarInit.derd[1].u)) * freeMotionScalarInit.angle_d_1 + (-sin(freeMotionScalarInit.derd[3].u)) * (freeMotionScalarInit.initAngle.R_rel.T[3,1] * cos(freeMotionScalarInit.derd[1].u) * freeMotionScalarInit.angle_d_1 + $DER.freeMotionScalarInit.initAngle.R_rel.T[3,1] * sin(freeMotionScalarInit.derd[1].u)) - cos(freeMotionScalarInit.derd[3].u) * freeMotionScalarInit.angle_d_3 * freeMotionScalarInit.initAngle.R_rel.T[3,1] * sin(freeMotionScalarInit.derd[1].u) - sin(freeMotionScalarInit.derd[3].u) * freeMotionScalarInit.angle_d_3 * cos(freeMotionScalarInit.derd[1].u)   [dynamic |0|0|0|0|]
19/19 (1): $DER.freeMotionScalarInit.initAngle.R_rel.T[1,3] = sin(freeMotionScalarInit.derd[3].u) * cos(freeMotionScalarInit.derd[1].u) * freeMotionScalarInit.angle_d_1 + cos(freeMotionScalarInit.derd[3].u) * freeMotionScalarInit.angle_d_3 * sin(freeMotionScalarInit.derd[1].u) + sin(freeMotionScalarInit.derd[3].u) * freeMotionScalarInit.angle_d_3 * freeMotionScalarInit.initAngle.R_rel.T[3,1] * cos(freeMotionScalarInit.derd[1].u) - cos(freeMotionScalarInit.derd[3].u) * (freeMotionScalarInit.initAngle.R_rel.T[3,1] * (-sin(freeMotionScalarInit.derd[1].u)) * freeMotionScalarInit.angle_d_1 + $DER.freeMotionScalarInit.initAngle.R_rel.T[3,1] * cos(freeMotionScalarInit.derd[1].u))   [dynamic |0|0|0|0|]
20/20 (1): $DER.freeMotionScalarInit.initAngle.R_rel.T[1,2] = sin(freeMotionScalarInit.derd[3].u) * (-sin(freeMotionScalarInit.derd[1].u)) * freeMotionScalarInit.angle_d_1 + cos(freeMotionScalarInit.derd[3].u) * freeMotionScalarInit.angle_d_3 * cos(freeMotionScalarInit.derd[1].u) + cos(freeMotionScalarInit.derd[3].u) * (freeMotionScalarInit.initAngle.R_rel.T[3,1] * cos(freeMotionScalarInit.derd[1].u) * freeMotionScalarInit.angle_d_1 + $DER.freeMotionScalarInit.initAngle.R_rel.T[3,1] * sin(freeMotionScalarInit.derd[1].u)) - sin(freeMotionScalarInit.derd[3].u) * freeMotionScalarInit.angle_d_3 * freeMotionScalarInit.initAngle.R_rel.T[3,1] * sin(freeMotionScalarInit.derd[1].u)   [dynamic |0|0|0|0|]
22/22 (1): freeMotionScalarInit.initAngle.R_rel.w[3] = freeMotionScalarInit.angle_d_3 + freeMotionScalarInit.initAngle.R_rel.T[3,1] * freeMotionScalarInit.angle_d_1   [dynamic |0|0|0|0|]

Equations where freeMotionScalarInit.angle_d_2 appears:

15/15 (1): $DER.freeMotionScalarInit.initAngle.R_rel.T[3,1] = cos(freeMotionScalarInit.derd[2].u) * freeMotionScalarInit.angle_d_2   [dynamic |0|0|0|0|]
18/18 (1): $DER.freeMotionScalarInit.initAngle.R_rel.T[2,1] = (-cos(freeMotionScalarInit.derd[3].u)) * freeMotionScalarInit.angle_d_3 * cos(freeMotionScalarInit.derd[2].u) - sin(freeMotionScalarInit.derd[3].u) * (-sin(freeMotionScalarInit.derd[2].u)) * freeMotionScalarInit.angle_d_2   [dynamic |0|0|0|0|]
21/21 (1): $DER.freeMotionScalarInit.initAngle.R_rel.T[1,1] = cos(freeMotionScalarInit.derd[3].u) * (-sin(freeMotionScalarInit.derd[2].u)) * freeMotionScalarInit.angle_d_2 - sin(freeMotionScalarInit.derd[3].u) * freeMotionScalarInit.angle_d_3 * cos(freeMotionScalarInit.derd[2].u)   [dynamic |0|0|0|0|]

Equations where both appear:

13/13 (1): $DER.freeMotionScalarInit.initAngle.R_rel.T[3,3] = cos(freeMotionScalarInit.derd[2].u) * (-sin(freeMotionScalarInit.derd[1].u)) * freeMotionScalarInit.angle_d_1 - sin(freeMotionScalarInit.derd[2].u) * freeMotionScalarInit.angle_d_2 * cos(freeMotionScalarInit.derd[1].u)   [dynamic |0|0|0|0|]
14/14 (1): $DER.freeMotionScalarInit.initAngle.R_rel.T[3,2] = sin(freeMotionScalarInit.derd[2].u) * freeMotionScalarInit.angle_d_2 * sin(freeMotionScalarInit.derd[1].u) - cos(freeMotionScalarInit.derd[2].u) * cos(freeMotionScalarInit.derd[1].u) * freeMotionScalarInit.angle_d_1   [dynamic |0|0|0|0|]
23/23 (1): freeMotionScalarInit.initAngle.R_rel.w[2] = cos(freeMotionScalarInit.derd[3].u) * freeMotionScalarInit.angle_d_2 + freeMotionScalarInit.initAngle.R_rel.T[2,1] * freeMotionScalarInit.angle_d_1   [dynamic |0|0|0|0|]
24/24 (1): freeMotionScalarInit.initAngle.R_rel.w[1] = sin(freeMotionScalarInit.derd[3].u) * freeMotionScalarInit.angle_d_2 + freeMotionScalarInit.initAngle.R_rel.T[1,1] * freeMotionScalarInit.angle_d_1   [dynamic |0|0|0|0|]

@Francesco could you tell me where the user expects these values to get solved when tagging them with StateSelect.never? That would help a lot! My default implementation tries to solve freeMotionScalarInit.angle_d_1 in equation 14 and freeMotionScalarInit.angle_d_2 in equation 21.

[casella]

Replying to Karim.Abdelhak:

1. The following models only break in MSL 3.1 and work in the current one. It seems like some min - max constraints don't work out the way they should. When running MixtureGases i realized that medium2.T equals medium2.T_degC which seems wrong, since T should be Kelvin and T_degC should be degrees Celsius. Do you know what exactly changed here from 3.1 to now?

I could double-check, but as far as I am concerned, MSL versions prior to 3.2.1 were affected by a lot of fundamental flaws, mainly regarding the fact that they were not really consistent with the standard. This is also the reason why trying to patch MSL 3.1 to solve these issues would be a lost cause.

Also, there should be no major issues upgrading existing models from MSL 3.x to MSL 3.2.3, so I won't be worried by this issue at all.

@Francesco Could you give me some thoughts on what i reported, especially to 1.? And is it possible that the selection of StateSelect.never states is actually correct or is it always wrong since there was a selection we could make without choosing them?

I'm not sure what you mean here. A variable with StateSelect.never should never be selected as a state. If this causes any problem, the problem should be reported and the translation aborted.

[casella]

Replying to Karim.Abdelhak:

@Francesco could you tell me where the user expects these values to get solved when tagging them with StateSelect.never? That would help a lot! My default implementation tries to solve freeMotionScalarInit.angle_d_1 in equation 14 and freeMotionScalarInit.angle_d_2 in equation 21.

I ran the model in Dymola, which was used to develop the MultiBody library. In the dsmodel.mof file I found

  freeMotionScalarInit.angle_d_2 := der(freeMotionScalarInit.initAngle.angle[2]);

// Dynamic selection of states.
  // State candidates
  // of which 1 shall be selected as a state,
  // and 7 shall not be a state:
  //   der(freeMotionScalarInit.initAngle.angle[2])
  //   der(freeMotionScalarInit.initAngle.angle[1])
  //   freeMotionScalarInit.initAngle.R_rel.w[3]
  //   freeMotionScalarInit.initAngle.R_rel.w[2]
  //   freeMotionScalarInit.initAngle.R_rel.w[1]
  //   bodyOfConstraint.body.w_a[3]
  //   bodyOfConstraint.body.w_a[2]
  //   bodyOfConstraint.body.w_a[1]
// The constraints form a linear system of equations, J*x = b,
// with J and b as given below.
  J[1, 1] := (joint.frame_a.R.T[3, 1]*constraint.frame_a.R.T[2, 1]+
    joint.frame_a.R.T[3, 2]*constraint.frame_a.R.T[2, 2]+joint.frame_a.R.T[3, 3]
    *constraint.frame_a.R.T[2, 3])*(cos(freeMotionScalarInit.initAngle.angle[1])
    *sin(freeMotionScalarInit.initAngle.angle[2]))-((joint.frame_a.R.T[1, 1]*
    constraint.frame_a.R.T[2, 1]+joint.frame_a.R.T[1, 2]*constraint.frame_a.R.T[2, 2]
    +joint.frame_a.R.T[1, 3]*constraint.frame_a.R.T[2, 3])*cos(freeMotionScalarInit.initAngle.angle[2])
    +(joint.frame_a.R.T[2, 1]*constraint.frame_a.R.T[2, 1]+joint.frame_a.R.T[2, 2]
    *constraint.frame_a.R.T[2, 2]+joint.frame_a.R.T[2, 3]*constraint.frame_a.R.T[2, 3])
    *(sin(freeMotionScalarInit.initAngle.angle[1])*sin(freeMotionScalarInit.initAngle.angle[2])));
  J[1, 2] := (joint.frame_a.R.T[2, 1]*constraint.frame_a.R.T[2, 1]+
    joint.frame_a.R.T[2, 2]*constraint.frame_a.R.T[2, 2]+joint.frame_a.R.T[2, 3]
    *constraint.frame_a.R.T[2, 3])*(cos(freeMotionScalarInit.initAngle.angle[2])
    *cos(freeMotionScalarInit.initAngle.angle[1]))+(joint.frame_a.R.T[3, 1]*
    constraint.frame_a.R.T[2, 1]+joint.frame_a.R.T[3, 2]*constraint.frame_a.R.T[2, 2]
    +joint.frame_a.R.T[3, 3]*constraint.frame_a.R.T[2, 3])*(cos(freeMotionScalarInit.initAngle.angle[2])
    *sin(freeMotionScalarInit.initAngle.angle[1]));
  J[2, 3] := -1.0;
  J[2, 6] := 1.0;
  J[3, 4] := -1.0;
  J[3, 7] := 1.0;
  J[4, 1] :=  -cos(freeMotionScalarInit.initAngle.angle[3]);
  J[4, 2] := sin(freeMotionScalarInit.initAngle.angle[3])*cos(freeMotionScalarInit.initAngle.angle[2]);
  J[4, 4] := 1.0;
  J[5, 5] := -1.0;
  J[5, 8] := 1.0;
  J[6, 1] :=  -sin(freeMotionScalarInit.initAngle.angle[3]);
  J[6, 2] :=  -cos(freeMotionScalarInit.initAngle.angle[3])*cos(freeMotionScalarInit.initAngle.angle[2]);
  J[6, 5] := 1.0;
  J[7, 1] := (joint.frame_a.R.T[1, 1]*constraint.frame_a.R.T[2, 1]+
    joint.frame_a.R.T[1, 2]*constraint.frame_a.R.T[2, 2]+joint.frame_a.R.T[1, 3]
    *constraint.frame_a.R.T[2, 3])*(cos(freeMotionScalarInit.initAngle.angle[3])
    *sin(freeMotionScalarInit.initAngle.angle[2]))-(joint.frame_a.R.T[2, 1]*
    constraint.frame_a.R.T[2, 1]+joint.frame_a.R.T[2, 2]*constraint.frame_a.R.T[2, 2]
    +joint.frame_a.R.T[2, 3]*constraint.frame_a.R.T[2, 3])*(sin(freeMotionScalarInit.initAngle.angle[1])
    *(cos(freeMotionScalarInit.initAngle.angle[3])*cos(freeMotionScalarInit.initAngle.angle[2])))
    +(joint.frame_a.R.T[3, 1]*constraint.frame_a.R.T[2, 1]+joint.frame_a.R.T[3, 2]
    *constraint.frame_a.R.T[2, 2]+joint.frame_a.R.T[3, 3]*constraint.frame_a.R.T[2, 3])
    *(cos(freeMotionScalarInit.initAngle.angle[1])*(cos(freeMotionScalarInit.initAngle.angle[3])
    *cos(freeMotionScalarInit.initAngle.angle[2])));
  J[7, 2] :=  -((cos(freeMotionScalarInit.initAngle.angle[2])*sin(
    freeMotionScalarInit.initAngle.angle[3])*(joint.frame_a.R.T[1, 1]*
    constraint.frame_a.R.T[2, 1]+joint.frame_a.R.T[1, 2]*constraint.frame_a.R.T[2, 2]
    +joint.frame_a.R.T[1, 3]*constraint.frame_a.R.T[2, 3])-(cos(freeMotionScalarInit.initAngle.angle[1])
    *cos(freeMotionScalarInit.initAngle.angle[3])-sin(freeMotionScalarInit.initAngle.angle[1])
    *(sin(freeMotionScalarInit.initAngle.angle[2])*sin(freeMotionScalarInit.initAngle.angle[3])))
    *(joint.frame_a.R.T[2, 1]*constraint.frame_a.R.T[2, 1]+joint.frame_a.R.T[2, 2]
    *constraint.frame_a.R.T[2, 2]+joint.frame_a.R.T[2, 3]*constraint.frame_a.R.T[2, 3])
    -(sin(freeMotionScalarInit.initAngle.angle[1])*cos(freeMotionScalarInit.initAngle.angle[3])
    +cos(freeMotionScalarInit.initAngle.angle[1])*(sin(freeMotionScalarInit.initAngle.angle[2])
    *sin(freeMotionScalarInit.initAngle.angle[3])))*(joint.frame_a.R.T[3, 1]*
    constraint.frame_a.R.T[2, 1]+joint.frame_a.R.T[3, 2]*constraint.frame_a.R.T[2, 2]
    +joint.frame_a.R.T[3, 3]*constraint.frame_a.R.T[2, 3]))*sin(freeMotionScalarInit.initAngle.angle[2])
    +(joint.frame_a.R.T[2, 1]*constraint.frame_a.R.T[2, 1]+joint.frame_a.R.T[2, 2]
    *constraint.frame_a.R.T[2, 2]+joint.frame_a.R.T[2, 3]*constraint.frame_a.R.T[2, 3])
    *(cos(freeMotionScalarInit.initAngle.angle[3])*sin(freeMotionScalarInit.initAngle.angle[2])
    *cos(freeMotionScalarInit.initAngle.angle[1])-sin(freeMotionScalarInit.initAngle.angle[3])
    *sin(freeMotionScalarInit.initAngle.angle[1]))+(joint.frame_a.R.T[3, 1]*
    constraint.frame_a.R.T[2, 1]+joint.frame_a.R.T[3, 2]*constraint.frame_a.R.T[2, 2]
    +joint.frame_a.R.T[3, 3]*constraint.frame_a.R.T[2, 3])*(sin(freeMotionScalarInit.initAngle.angle[3])
    *cos(freeMotionScalarInit.initAngle.angle[1])+cos(freeMotionScalarInit.initAngle.angle[3])
    *sin(freeMotionScalarInit.initAngle.angle[2])*sin(freeMotionScalarInit.initAngle.angle[1])));
  J[7, 3] := cos(freeMotionScalarInit.initAngle.angle[2])*sin(freeMotionScalarInit.initAngle.angle[3])
    *(joint.frame_a.R.T[1, 1]*constraint.frame_a.R.T[2, 1]+joint.frame_a.R.T[1, 2]
    *constraint.frame_a.R.T[2, 2]+joint.frame_a.R.T[1, 3]*constraint.frame_a.R.T[2, 3])
    -(cos(freeMotionScalarInit.initAngle.angle[1])*cos(freeMotionScalarInit.initAngle.angle[3])
    -sin(freeMotionScalarInit.initAngle.angle[1])*(sin(freeMotionScalarInit.initAngle.angle[2])
    *sin(freeMotionScalarInit.initAngle.angle[3])))*(joint.frame_a.R.T[2, 1]*
    constraint.frame_a.R.T[2, 1]+joint.frame_a.R.T[2, 2]*constraint.frame_a.R.T[2, 2]
    +joint.frame_a.R.T[2, 3]*constraint.frame_a.R.T[2, 3])-(sin(freeMotionScalarInit.initAngle.angle[1])
    *cos(freeMotionScalarInit.initAngle.angle[3])+cos(freeMotionScalarInit.initAngle.angle[1])
    *(sin(freeMotionScalarInit.initAngle.angle[2])*sin(freeMotionScalarInit.initAngle.angle[3])))
    *(joint.frame_a.R.T[3, 1]*constraint.frame_a.R.T[2, 1]+joint.frame_a.R.T[3, 2]
    *constraint.frame_a.R.T[2, 2]+joint.frame_a.R.T[3, 3]*constraint.frame_a.R.T[2, 3]);
  b[2] := freeMotionScalarInit.initAngle.R_rel.T[3, 1]*joint.frame_a.R.w[1]+
    freeMotionScalarInit.initAngle.R_rel.T[3, 2]*joint.frame_a.R.w[2]+
    freeMotionScalarInit.initAngle.R_rel.T[3, 3]*joint.frame_a.R.w[3];
  b[3] := freeMotionScalarInit.initAngle.R_rel.T[2, 1]*joint.frame_a.R.w[1]+
    freeMotionScalarInit.initAngle.R_rel.T[2, 2]*joint.frame_a.R.w[2]+
    freeMotionScalarInit.initAngle.R_rel.T[2, 3]*joint.frame_a.R.w[3];
  b[5] := freeMotionScalarInit.initAngle.R_rel.T[1, 1]*joint.frame_a.R.w[1]+
    freeMotionScalarInit.initAngle.R_rel.T[1, 2]*joint.frame_a.R.w[2]+
    freeMotionScalarInit.initAngle.R_rel.T[1, 3]*joint.frame_a.R.w[3];
algorithm // Torn part
  der(freeMotionScalarInit.initAngle.R_rel.T[3, 1]) := der(freeMotionScalarInit.initAngle.angle[2])
    *cos(freeMotionScalarInit.initAngle.angle[2]);
  der(freeMotionScalarInit.initAngle.R_rel.T[3, 2]) := der(freeMotionScalarInit.initAngle.angle[2])
    *sin(freeMotionScalarInit.initAngle.angle[2])*sin(freeMotionScalarInit.initAngle.angle[1])
    -cos(freeMotionScalarInit.initAngle.angle[2])*(der(freeMotionScalarInit.initAngle.angle[1])
    *cos(freeMotionScalarInit.initAngle.angle[1]));
  der(freeMotionScalarInit.initAngle.R_rel.T[3, 3]) :=  -(der(freeMotionScalarInit.initAngle.angle[2])
    *sin(freeMotionScalarInit.initAngle.angle[2])*cos(freeMotionScalarInit.initAngle.angle[1])
    +cos(freeMotionScalarInit.initAngle.angle[2])*(der(freeMotionScalarInit.initAngle.angle[1])
    *sin(freeMotionScalarInit.initAngle.angle[1])));
  der(constraint.frame_b.R.T[3, 1]) := der(freeMotionScalarInit.initAngle.R_rel.T[3, 1])
    *joint.frame_a.R.T[1, 1]+der(freeMotionScalarInit.initAngle.R_rel.T[3, 2])*
    joint.frame_a.R.T[2, 1]+der(freeMotionScalarInit.initAngle.R_rel.T[3, 3])*
    joint.frame_a.R.T[3, 1];
  der(constraint.frame_b.R.T[3, 2]) := der(freeMotionScalarInit.initAngle.R_rel.T[3, 1])
    *joint.frame_a.R.T[1, 2]+der(freeMotionScalarInit.initAngle.R_rel.T[3, 2])*
    joint.frame_a.R.T[2, 2]+der(freeMotionScalarInit.initAngle.R_rel.T[3, 3])*
    joint.frame_a.R.T[3, 2];
  der(constraint.frame_b.R.T[3, 3]) := der(freeMotionScalarInit.initAngle.R_rel.T[3, 1])
    *joint.frame_a.R.T[1, 3]+der(freeMotionScalarInit.initAngle.R_rel.T[3, 2])*
    joint.frame_a.R.T[2, 3]+der(freeMotionScalarInit.initAngle.R_rel.T[3, 3])*
    joint.frame_a.R.T[3, 3];
  der(freeMotionScalarInit.initAngle.angle[3]) := freeMotionScalarInit.initAngle.R_rel.w[3]
    -sin(freeMotionScalarInit.initAngle.angle[2])*der(freeMotionScalarInit.initAngle.angle[1]);
  der(freeMotionScalarInit.initAngle.R_rel.T[1, 1]) :=  -(der(freeMotionScalarInit.initAngle.angle[3])
    *sin(freeMotionScalarInit.initAngle.angle[3])*cos(freeMotionScalarInit.initAngle.angle[2])
    +cos(freeMotionScalarInit.initAngle.angle[3])*(der(freeMotionScalarInit.initAngle.angle[2])
    *sin(freeMotionScalarInit.initAngle.angle[2])));
  der(freeMotionScalarInit.initAngle.R_rel.T[1, 2]) := der(freeMotionScalarInit.initAngle.angle[3])
    *cos(freeMotionScalarInit.initAngle.angle[3])*cos(freeMotionScalarInit.initAngle.angle[1])
    -sin(freeMotionScalarInit.initAngle.angle[3])*(der(freeMotionScalarInit.initAngle.angle[1])
    *sin(freeMotionScalarInit.initAngle.angle[1]))+(cos(freeMotionScalarInit.initAngle.angle[3])
    *(der(freeMotionScalarInit.initAngle.angle[2])*cos(freeMotionScalarInit.initAngle.angle[2]))
    -der(freeMotionScalarInit.initAngle.angle[3])*sin(freeMotionScalarInit.initAngle.angle[3])
    *sin(freeMotionScalarInit.initAngle.angle[2]))*sin(freeMotionScalarInit.initAngle.angle[1])
    +cos(freeMotionScalarInit.initAngle.angle[3])*sin(freeMotionScalarInit.initAngle.angle[2])
    *(der(freeMotionScalarInit.initAngle.angle[1])*cos(freeMotionScalarInit.initAngle.angle[1]));
  der(freeMotionScalarInit.initAngle.R_rel.T[1, 3]) := der(freeMotionScalarInit.initAngle.angle[3])
    *cos(freeMotionScalarInit.initAngle.angle[3])*sin(freeMotionScalarInit.initAngle.angle[1])
    +sin(freeMotionScalarInit.initAngle.angle[3])*(der(freeMotionScalarInit.initAngle.angle[1])
    *cos(freeMotionScalarInit.initAngle.angle[1]))-(cos(freeMotionScalarInit.initAngle.angle[3])
    *(der(freeMotionScalarInit.initAngle.angle[2])*cos(freeMotionScalarInit.initAngle.angle[2]))
    -der(freeMotionScalarInit.initAngle.angle[3])*sin(freeMotionScalarInit.initAngle.angle[3])
    *sin(freeMotionScalarInit.initAngle.angle[2]))*cos(freeMotionScalarInit.initAngle.angle[1])
    +cos(freeMotionScalarInit.initAngle.angle[3])*sin(freeMotionScalarInit.initAngle.angle[2])
    *(der(freeMotionScalarInit.initAngle.angle[1])*sin(freeMotionScalarInit.initAngle.angle[1]));
  der(constraint.frame_b.R.T[1, 1]) := der(freeMotionScalarInit.initAngle.R_rel.T[1, 1])
    *joint.frame_a.R.T[1, 1]+der(freeMotionScalarInit.initAngle.R_rel.T[1, 2])*
    joint.frame_a.R.T[2, 1]+der(freeMotionScalarInit.initAngle.R_rel.T[1, 3])*
    joint.frame_a.R.T[3, 1];
  der(constraint.frame_b.R.T[1, 2]) := der(freeMotionScalarInit.initAngle.R_rel.T[1, 1])
    *joint.frame_a.R.T[1, 2]+der(freeMotionScalarInit.initAngle.R_rel.T[1, 2])*
    joint.frame_a.R.T[2, 2]+der(freeMotionScalarInit.initAngle.R_rel.T[1, 3])*
    joint.frame_a.R.T[3, 2];
  der(constraint.frame_b.R.T[1, 3]) := der(freeMotionScalarInit.initAngle.R_rel.T[1, 1])
    *joint.frame_a.R.T[1, 3]+der(freeMotionScalarInit.initAngle.R_rel.T[1, 2])*
    joint.frame_a.R.T[2, 3]+der(freeMotionScalarInit.initAngle.R_rel.T[1, 3])*
    joint.frame_a.R.T[3, 3];
// End dynamic selection of states

Does this help you?

[Karim.Abdelhak]

Replying to casella:

Does this help you?

From what i can see here, those are actually part of the dynamic state selection set dispite beeing tagged with StateSelect.never. It seems like dymola sets up the set with the lower order derivatives der(freeMotionScalarInit.initAngle.angle[2]) as states instead of using the states, those derivatives don't inherit the StateSelect attribute of the 'alias'. We generally use a different approach here, but the information is present. So maybe we should just not force StateSelect.never on states which have a direct representation as the derivative of another state.

Seems like a workaround but should result in the same sets. I will try that tomorrow, thank you Francesco!

[Karim.Abdelhak]

It works out and i also understand now why it has to be that way. When simulating the states which have a direct representation as a derivative of another, we naturally would simulate the derivatives of the other, since the states are known. Even if the state itself can't be picked as a state to be simulated, the derivative of the other (which is an alias) still could be.

It is a bit convoluted, since we still simulate the state with StateSelect.never instead of the derivative of the other, but that's mostly due to the way we derive expressions. It is the same as dymola but not that intuitive. I shortly explained it in a inline comment in the code and linked it to this ticket.

The last remaining problems are some linearizations from optimica:

tests / testsuite-gcc / testSteamPipe.mos – openmodelica_linearization21s
tests / testsuite-gcc / testDrumBoiler.mos – openmodelica_linearization21s

The easy way would be to disable my changes for linearization optimization problems. The solver converges to a different linearization which should be wrong. It reports that the jacobian is nonlinear and that could lead to errors, i guess this case is just not implemented. I actually don't know that much about our optimization implementation, so that would be quite the task.

@Francesco Any opinion on this?

[casella]

Replying to Karim.Abdelhak:

From what i can see here, those are actually part of the dynamic state selection set dispite beeing tagged with StateSelect.never. It seems like dymola sets up the set with the lower order derivatives der(freeMotionScalarInit.initAngle.angle[2]) as states instead of using the states, those derivatives don't inherit the StateSelect attribute of the 'alias'.

I checked the Specification, sect. 4.8.8.1, there is no mention of 'alias variables', so I guess we should interpret it strictly.

[casella]

Replying to Karim.Abdelhak:

The last remaining problems are some linearizations from optimica:

tests / testsuite-gcc / testSteamPipe.mos – openmodelica_linearization21s
tests / testsuite-gcc / testDrumBoiler.mos – openmodelica_linearization21s

If you are referring to

​https://github.com/OpenModelica/OpenModelica-testsuite/blob/master/openmodelica/linearization/testDrumBoiler.mos

I can't see any reference to optimization. Is this the test you refer to?

Linearization per se has nothing to do with optimization, and is an essential feature that I personally use a lot to support control system design. It should definitely not break. However, as it is well-known, there are infinitely many equivalent linearizations of linear systems, e.g. by changing the order of the states in the state vector, which is arbitrary. What should be the same (up to rounding errors) is the transfer function between inputs and outputs. Did you check that?

[Karim.Abdelhak]

Replying to casella:

Replying to Karim.Abdelhak:

The last remaining problems are some linearizations from optimica:

tests / testsuite-gcc / testSteamPipe.mos – openmodelica_linearization21s
tests / testsuite-gcc / testDrumBoiler.mos – openmodelica_linearization21s

If you are referring to

​https://github.com/OpenModelica/OpenModelica-testsuite/blob/master/openmodelica/linearization/testDrumBoiler.mos

I can't see any reference to optimization. Is this the test you refer to?

Linearization per se has nothing to do with optimization, and is an essential feature that I personally use a lot to support control system design. It should definitely not break. However, as it is well-known, there are infinitely many equivalent linearizations of linear systems, e.g. by changing the order of the states in the state vector, which is arbitrary. What should be the same (up to rounding errors) is the transfer function between inputs and outputs. Did you check that?

This is the log:

==== Log /tmp/omc-rtest-unknown/openmodelica/linearization/testDrumBoiler.mos_temp9985/log-testDrumBoiler.mos
true
""
true
true
record SimulationResult
resultFile = "testDrumBoilerLin_res.mat",
simulationOptions = "startTime = 0.0, stopTime = 0.0, numberOfIntervals = 500, tolerance = 1e-06, method = 'dassl', fileNamePrefix = 'testDrumBoilerLin', options = '', outputFormat = 'mat', variableFilter = '.*', cflags = '', simflags = ''",
messages = "stdout | info | Linearization will performed at point of time: 0.000000
LOG_SUCCESS | info | The initialization finished successfully without homotopy method.
LOG_SUCCESS | info | The simulation finished successfully.
stdout | info | Linear model is created!
"
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.
Error: Jacobian A contains non-linear components. This indicates a singular system or internal generation errors.
"
true
"model linear_testDrumBoilerLin
parameter Integer n = 3 \"number of states\";
parameter Integer p = 2 \"number of inputs\";
parameter Integer q = 4 \"number of outputs\";
parameter Real x0[n] = {0, 82766887470.64922, 364248.1518687711};
parameter Real u0[p] = {0, 0};
parameter Real A[n, n] = [0, 0, -8.293549697411924e-06; 5000000, 0, 0; 10, 0, 0];
parameter Real B[n, p] = [0, 0; -2674949.557221949, 1000000; -1, 0];
parameter Real C[q, n] = [0, 0, -0.0005049700197164287; 0, 0, 0.000995225963689431; 0, 0, -3.525747112297686e-05; 0, 0, 0];
parameter Real D[q, p] = [0, 0; 0, 0; 0, 0; 1, 0];
Real x[n](start = x0);
input Real u[p](start = u0);
output Real y[q];
Real 'x_controller.x' = x[1];
Real 'x_evaporator.U' = x[2];
Real 'x_evaporator.m' = x[3];
Real 'u_Y_Valve' = u[1];
Real 'u_q_F' = u[2];
Real 'y_T_S' = y[1];
Real 'y_V_l' = y[2];
Real 'y_p_S' = y[3];
Real 'y_qm_S' = y[4];
equation
der(x) = A * x + B * u;
y = C * x + D * u;
end linear_testDrumBoilerLin;"
Equation mismatch: diff says:
--- /tmp/omc-rtest-unknown/openmodelica/linearization/testDrumBoiler.mos_temp9985/equations-expected2019-08-21 12:55:29.891416456 +0000
+++ /tmp/omc-rtest-unknown/openmodelica/linearization/testDrumBoiler.mos_temp9985/equations-got2019-08-21 12:55:50.431342292 +0000
@@ -10,28 +10,29 @@
LOG_SUCCESS | info | The simulation finished successfully.
stdout | info | Linear model is created!
"
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.
+Error: Jacobian A contains non-linear components. This indicates a singular system or internal generation errors.
"
true
"model linear_testDrumBoilerLin
parameter Integer n = 3 \"number of states\";
parameter Integer p = 2 \"number of inputs\";
parameter Integer q = 4 \"number of outputs\";
-parameter Real x0[n] = {0, 67, 100000};
+parameter Real x0[n] = {0, 82766887470.64922, 364248.1518687711};
parameter Real u0[p] = {0, 0};
-parameter Real A[n, n] = [-0, -0.008333333333333333, -0; 0.01043953430921842, -0.01043953430921842, 0; 0.1178989396709848, -0.1178989396709848, 4.135580766728708e-15];
-parameter Real B[n, p] = [0, 0; -0.001308242749261165, 0.0001170710024614632; -19.14622757506175, 8.475892309753029];
-parameter Real C[q, n] = [0, 0, 0.0001439468903880623; 0, 1, 0; 0, 0, 1e-05; 0, 0, 0];
+parameter Real A[n, n] = [0, 0, -8.293549697411924e-06; 5000000, 0, 0; 10, 0, 0];
+parameter Real B[n, p] = [0, 0; -2674949.557221949, 1000000; -1, 0];
+parameter Real C[q, n] = [0, 0, -0.0005049700197164287; 0, 0, 0.000995225963689431; 0, 0, -3.525747112297686e-05; 0, 0, 0];
parameter Real D[q, p] = [0, 0; 0, 0; 0, 0; 1, 0];
Real x[n](start = x0);
input Real u[p](start = u0);
output Real y[q];
Real 'x_controller.x' = x[1];
-Real 'x_evaporator.V_l' = x[2];
-Real 'x_evaporator.p' = x[3];
+Real 'x_evaporator.U' = x[2];
+Real 'x_evaporator.m' = x[3];
Real 'u_Y_Valve' = u[1];
Real 'u_q_F' = u[2];
Real 'y_T_S' = y[1];
Real 'y_V_l' = y[2];
Real 'y_p_S' = y[3];
Equation mismatch: omc-diff says:
Failed '"' 'E'
Line 15: Text differs:
expected: "
got: Error: Jacobian A contains non
== 1 out of 1 tests failed [openmodelica/linearization/testDrumBoiler.mos_temp9985, time: 21]

It does not seem correct.

[casella]

Replying to Karim.Abdelhak:

-Real 'x_evaporator.V_l' = x[2];
-Real 'x_evaporator.p' = x[3];
+Real 'x_evaporator.U' = x[2];
+Real 'x_evaporator.m' = x[3];

Something obviously went wrong with state selection. If you look at the ​source code of the evaporator model, U and m are the variables that appear as arguments of der() operators in the model, and they have StateSelect.default, while p and V_l have StateSelect.prefer. This should be enough for the latter to be chosen as states, which is actually what happens if the model is simulated.

I'm not sure why this doesn't happen when linearizing the model, but this should be investigated. In principle, there is no difference between generating the code for linearization and generating the code for simulation with an ODE solver.

Do you have any idea why your PR could affect the linearize operation regarding state selection?

[Karim.Abdelhak]

Replying to casella:

Replying to Karim.Abdelhak:

-Real 'x_evaporator.V_l' = x[2];
-Real 'x_evaporator.p' = x[3];
+Real 'x_evaporator.U' = x[2];
+Real 'x_evaporator.m' = x[3];

Something obviously went wrong with state selection. If you look at the ​source code of the evaporator model, U and m are the variables that appear as arguments of der() operators in the model, and they have StateSelect.default, while p and V_l have StateSelect.prefer. This should be enough for the latter to be chosen as states, which is actually what happens if the model is simulated.

I'm not sure why this doesn't happen when linearizing the model, but this should be investigated. In principle, there is no difference between generating the code for linearization and generating the code for simulation with an ODE solver.

Do you have any idea why your PR could affect the linearize operation regarding state selection?

Yes, i think i have quite the idea now. So what i basically changed is that only variables which appear derived in the system or those which have StateSelect.always are considered states. Before also variables with StateSelect.prefer were considered states, even if they don't appear derived. In many cases this previously forced state selection on the system, where none should be needed. The variables tagged with StateSelect.prefer did not appear derived anywhere, so there was no equation for the derivative to be solved in. From deriving all equations where these variables appeared (and other dependent equations) artificial state sets were formed which allowed those variables to become chosen as states.

In principle i understand why this was implemented in this way, but from my understanding a derivative of a variable that does not appear naturally should only be forced in the system when tagged with StateSelect.always. In the case of our heating system we had the problem, that some derivatives with StateSelect.prefer were forced in and created these artificial state sets with equations that should not be in there.

Now that i changed it, some of the variables which have StateSelect.prefer and don't appear derived are just considered algebraic variables. I don't see another way of making systems like the heating system work, where people just tag stuff with StateSelect attributes without knowing if they are actually states.

I think all StateSelect attributes besides StateSelect.always should not have any impact on non-state variables or systems that don't even have state selection naturally.

Which is the way to go?

[casella]

@Karim,

the way I always understood it is that StateSelect.always is a kind of doomsday device, like Mr. Draghi's bazooka (remember the "whatever it takes" quote?). It may be used, and if fact it is used in the MultiBody library, but it can easily lead to situations with more StateSelect.always variables than states after index reduction, which prevent continuing with the code generation. In fact, you see this situation explicitly mentioned, e.g., in the ​Rotational joint model.

This situation is always possible when assembling object-oriented models.

For Fluid models, the preferred states are usually pressure and temperature or enthalpy - however, it is often the case that pressure constraints are generated by connecting fluid ports, so using StateSelect.always on pressures is in general a very bad idea, because it may lead to compilation failures in those (otherwise perfectly legitimate) cases. This situation should be handled automatically, without the need by the modeller to expliticly disable the StateSelect.always attribute in such cases.

In normal situations, one should use StateSelect.prefer to specify that he/she wants a certain variable to become a state. As I see it, reading the specification, whether or not it appears diffentiated, it shouldn't matter at all. All people doing fluid modelling in Modelica expect variables with StateSelect.prefer that are not differentiated to become states. In fact, there is even a flag preferredMediumStates in the basic definition of the ​BasicProperties model for this purpose. It would be weird if OMC ignores it.

The tool should then try to make the StateSelect.prefer variables become states. Of course, in case there are more StateSelect.prefer variables than states after index reduction, see, e.g., the HeatingSystem case earlier in this ticket, then it is perfectly fine that some of them are not selected as states. In fact, that's the whole point of using prefer instead of always, which is too rigid. In this case, the tool should pick the states that lead to the easiest-to-solve model structure. But I don't think it's fine that we end up having a variable with StateSelect.prefer which is not a state and another variable with StateSelect.default that is. At least, this is not what I expected since I started doing fluid modelling in Modelica circa 2003 :)

It may even possible that selecting a StateSelect.prefer variable as state requires to solve implicit equations to solve for the original differentiated variables (this is not the case here). Even in this case, making this selection may have other beneficial effects, e.g. avoiding the solution of difficult implicit systems in some other parts of the model, so I would tend to favour the modeller's knowledge embedded in the StateSelect attribute and go for it anyway. But this is my interpretation and may be questionable.

Summing up, my recommendation for StateSelect.prefer variables is that the tool tries hard enough (to which extent is a matter of debate) to make them states, whether or not they are differentiated in the original model. In case there are more StateSelect.prefer variable than states in the index-1 reduced model, then the tool is free to select a subset that leads to the most efficient code at runtime.

[casella]

In any case, I still wonder why we get a different state selection for simulation and for linearization - shouldn't it be the same?

[Karim.Abdelhak]

Replying to casella:

In any case, I still wonder why we get a different state selection for simulation and for linearization - shouldn't it be the same?

Yes indeed, the one i provided is the one where i prevented the StateSelect.prefer variables from being forced as states. They are the same for linearization and simulation.

I guess i have to dig deeper into this problem once again since there seems to be something else. As you stated in this ticket, dymola reports

We cannot differentiate equations for the following variables
having stateSelect = StateSelect.prefer.
They will be treated as if they had stateSelect = StateSelect.default:
pump.medium.h
pump.medium.p

I don't know how we manage to derive these equations and dymola does not, maybe this 'deranking' of StateSelect attributes is what we need.

[Karim.Abdelhak]

I implemented some analytics and tried to reproduce what dymola is reporting. I got to the following subset of constraint equations which can't be derived by all candidates and the candidates in question:

 candidates not usable for partial differentiation   (8)
 ========================================
 tank.medium.h 
 heater.mediums[1].h
 heater.mediums[1].p
 pump.medium.p
 pump.medium.h
 radiator.mediums[1].h
 pipe.mediums[1].h
 pipe.mediums[2].h
 
 not partially differentiable equations  (8)
 ========================================
 1/1 (1): tank.portInDensities[2] = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(tank.vessel_ps_static[2], tank.medium.h, 0, 0).rho   [dynamic |0|0|0|0|]
 2/2 (1): heater.heatTransfer.Ts[1] = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(heater.mediums[1].p, heater.mediums[1].h, heater.statesFM[2].phase, 0).T   [dynamic |0|0|0|0|]
 3/3 (1): heater.flowModel.mus[2] = Modelica.Media.Water.IF97_Utilities.dynamicViscosity(heater.statesFM[2].d, heater.heatTransfer.Ts[1], heater.mediums[1].p, heater.statesFM[2].phase)   [dynamic |0|0|0|0|]
 4/4 (1): heater.statesFM[1].T = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(pump.medium.p, pump.medium.h, 0, 0).T   [dynamic |0|0|0|0|]
 5/5 (1): heater.statesFM[1].d = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(pump.medium.p, pump.medium.h, 0, 0).rho   [dynamic |0|0|0|0|]
 6/6 (1): radiator.statesFM[2].d = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(radiator.mediums[1].p, radiator.mediums[1].h, radiator.statesFM[2].phase, 0).rho   [dynamic |0|0|0|0|]
 7/7 (1): pipe.statesFM[1].d = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(pipe.mediums[1].p, pipe.mediums[1].h, pipe.statesFM[1].phase, 0).rho   [dynamic |0|0|0|0|]
 8/8 (1): pipe.statesFM[2].d = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(pipe.mediums[2].p, pipe.mediums[2].h, pipe.statesFM[2].phase, 0).rho   [dynamic |0|0|0|0|]

It seems like only the function calls from Modelica.Media.Water.IF97_Utilities cause problems, which seems reasonable.

I don't know why dymola only detects

 pump.medium.p
 pump.medium.h

as states with StateSelect.prefer for which we cannot differentiate. It looks like equation 4 and 5 are the ones which produce this.

I could now try to revoke the decision on these variables to be states and make them algebraic variables again, that would result in changing the original system and having to match it again. Most probably i can salvage the previous matching, but it's still very costly to do that.

With this approach it would most probably work, but still differ from dymolas result since i am revoking more state selection preferences. Also the static choice of states would differ, some states which are chosen by dymola would be algebraic variables in my approach. Most probably following equation would have the same problems and can't be differentiated by all candidates, once i have fixed if-equations:

0.0 = if tank.regularFlow[2] then tank.ports[2].p - homotopy(tank.vessel_ps_static[2] + 0.5 * tank.portAreas[2] ^ (-2.0) * Modelica.Fluid.Utilities.regSquare2(-m_flow, tank.m_flow_small, (-1.0 + tank.portsData[2].zeta_in  + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2], (1.0 +      tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.ports_penetration[2] * tank.     heatTransfer.states[1].d), false, 1.0), tank.vessel_ps_static[2]) else if tank.inFlow[2] then tank.ports[2].p - tank.vessel_ps_static[2] else -m_flow

We talked about this one and that it should not be differentiated.

Since this is much work i would like to hear your opinion @Francesco. My basic problem is that i don't know on which information the decision from dymola is made which StateSelect.prefer attributes should be ignored and which should not. From what i can see here there is no difference between any of them.

Info: All of the shown equations can be differentiated by another state candidate, only the if equation can't be matched at all.

[casella]

Replying to Karim.Abdelhak:

I implemented some analytics and tried to reproduce what dymola is reporting. I got to the following subset of constraint equations which can't be derived by all candidates and the candidates in question:

not partially differentiable equations (8)
========================================
1/1 (1): tank.portInDensities[2] = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(tank.vessel_ps_static[2], tank.medium.h, 0, 0).rho [dynamic |0|0|0|0|]
2/2 (1): heater.heatTransfer.Ts[1] = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(heater.mediums[1].p, heater.mediums[1].h, heater.statesFM[2].phase, 0).T [dynamic |0|0|0|0|]
3/3 (1): heater.flowModel.mus[2] = Modelica.Media.Water.IF97_Utilities.dynamicViscosity(heater.statesFM[2].d, heater.heatTransfer.Ts[1], heater.mediums[1].p, heater.statesFM[2].phase) [dynamic |0|0|0|0|]
4/4 (1): heater.statesFM[1].T = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(pump.medium.p, pump.medium.h, 0, 0).T [dynamic |0|0|0|0|]
5/5 (1): heater.statesFM[1].d = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(pump.medium.p, pump.medium.h, 0, 0).rho [dynamic |0|0|0|0|]
6/6 (1): radiator.statesFM[2].d = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(radiator.mediums[1].p, radiator.mediums[1].h, radiator.statesFM[2].phase, 0).rho [dynamic |0|0|0|0|]
7/7 (1): pipe.statesFM[1].d = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(pipe.mediums[1].p, pipe.mediums[1].h, pipe.statesFM[1].phase, 0).rho [dynamic |0|0|0|0|]
8/8 (1): pipe.statesFM[2].d = Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(pipe.mediums[2].p, pipe.mediums[2].h, pipe.statesFM[2].phase, 0).rho [dynamic |0|0|0|0|]

It seems like only the function calls from Modelica.Media.Water.IF97_Utilities cause problems, which seems reasonable.

I guess that's where the problem is. In order to change the states from m,U to p,H, it is necessary to get the derivatives of density w.r.t. p and h. The Modelica.Media IF97 model has a fairly intricate mechanism bases on Inline/LateInline annotations in order to enable a sort of caching mechanism, to avoid repeated calls to the core state equations when computing different properties, while at the same time providing derivative information.

Take, for example, the tank model, which is derived from Modelica.Fluid.Vessels.BaseClasses.PartialLumpedVessel. In there, you find this equation

    portInDensities[i] = Medium.density(Medium.setState_phX(vessel_ps_static[i], inStream(ports[i].h_outflow), inStream(ports[i].Xi_outflow)));

where Medium = Modelica.Media.Water.StandardWater = Modelica.Media.Water.WaterIF97_ph.

Function density returns state.d and is expected to be inlined.

state is the output of function setState_phX, which contains

  algorithm
    state := ThermodynamicState(
          d=density_ph(
            p,
            h,
            region=region),
          T=temperature_ph(
            p,
            h,
            region=region),
          phase=if region == 0 then 0 else if region==4 then 2 else 1,
          h=h,
          p=p);
    annotation (Inline=true);

hence portInDensities[i]=density_ph(...)

In turn density_ph() has this body:

  algorithm
    d := IF97_Utilities.rho_ph(
          p,
          h,
          phase,
          region);
    annotation (Inline=true);

so portInDensities[i]=IF97_Utilities.rho_ph(...), which has this body:

algorithm
  rho := rho_props_ph(
      p,
      h,
      waterBaseProp_ph(
        p,
        h,
        phase,
        region));
  annotation (Inline=true);

And therein lies the rub. The function rho_props_ph does the dirty trick:

algorithm
  rho := properties.rho;
  annotation (
    derivative(noDerivative=properties) = rho_ph_der,
    Inline=false,
    LateInline=true);

it returns the field d of the record properties, which is the output of the function waterBaseProp_ph, so it doesn't do any actual computation, but it has a derivative annotation, and a Inline = false, LateInline = true annotations.

Apparently, something went wrong there. When the index reduction analysis is performed, the call to rho_props_ph should not be inlined yet, as it instead is, since the call to waterBaseProp_ph shows up in the log. Hence, the derivative should be available thanks to the derivative annotation. The inlining of rho_props_ph should be done afterwards, once the index reduction step is complete.

I understand that is the reason why dymola doesn't (correctly) detect those variables as not differentiable.

I thought this machinery had been sorted out for good a few years ago, but it may be the case that there are still some wrinkles that need to be ironed out. @Karim, can you please check?

[Karim.Abdelhak]

Replying to casella:

Apparently, something went wrong there. When the index reduction analysis is performed, the call to rho_props_ph should not be inlined yet, as it instead is, since the call to waterBaseProp_ph shows up in the log. Hence, the derivative should be available thanks to the derivative annotation. The inlining of rho_props_ph should be done afterwards, once the index reduction step is complete.

I understand that is the reason why dymola doesn't (correctly) detect those variables as not differentiable.

I thought this machinery had been sorted out for good a few years ago, but it may be the case that there are still some wrinkles that need to be ironed out. @Karim, can you please check?

I found the origin of the problem and prevented premature inlining for this case, but the problem seems to remain. Not inlining the function preserves the information about differentiation w.r.t. time, but in this case we are trying to differentiate w.r.t. all state candidates. That also seems to be what dymola is doing since they report

We cannot differentiate equations for the following variables
having stateSelect = StateSelect.prefer.
They will be treated as if they had stateSelect = StateSelect.default:
pump.medium.h
pump.medium.p

We don't have any information on partial derivatives for those functions right?

Even if we had information on partial derivatives of functions, i don't see how the two state candidates dymola is reporting, differ from the others. Is there a way in dymola to report which equations could not be differentiated by those variables?

[casella]

Replying to Karim.Abdelhak:

I found the origin of the problem and prevented premature inlining for this case

Good.

but the problem seems to remain. Not inlining the function preserves the information about differentiation w.r.t. time, but in this case we are trying to differentiate w.r.t. all state candidates.

I beg your pardon, but this thread is becoming a bit long and I'm afraid I lost track of some bits and pieces. What is exactly the reason why you want to differentiate w.r.t. each variable individually?

What I understand is that we have a system like this

der(m) = aaa;
der(U) = bbb;
m = rho*V;
rho = rho_ph(p,h);
U = m*u;
u = h - p*V;

which has original two states m and U and four algebraic equations relating them to the variables with StateSelect.prefer, i.e., p and h.
if I differentiate those algebraic constraints, I get

der(m) = aaa;
der(U) = bbb;
der(m) = der(rho)*V;
der(rho) = rho_ph_der(p,h, der(p), der(h)); // thanks to the derivative annotation
der(U) = der(m)*u + m*der(u);
der(u) = der(h) - der(p)*V;

Now, to enforce StateSelect.prefer on p and h, you want der(p) and der(h) to show up for code generation. Hence, you have to mark all other der's in the set as dummies. That's it. I'm not sure how this maps on your algorithms using MSSS's, but why would you need individual partial derivatives of rho_ph w.r.t. each of its arguments? I can't see any use for them.

That also seems to be what dymola is doing since they report
We cannot differentiate equations for the following variables
having stateSelect = StateSelect.prefer.
They will be treated as if they had stateSelect = StateSelect.default:
pump.medium.h
pump.medium.p

This case is different. The pump model has no dynamic states, so there is no state variable change to be made here, because there are no states at all. Again, I'm not sure how this maps into your algorithms, but as I see it, that's where the difference is.

[Karim.Abdelhak]

Replying to casella:

I beg your pardon, but this thread is becoming a bit long and I'm afraid I lost track of some bits and pieces. What is exactly the reason why you want to differentiate w.r.t. each variable individually?

Yes this gets quite convoluted. The reason we want to differentiate w.r.t. the candidates is that we look for solvability. We need to ensure that the jacobian of the constraint equations / dummy states is regular at all times. Since we have more candidates than equations we try to chose subset of states to be chosen as dummys such that the resulting jacobian is constant and regular (static index reduction). If such a set does not exist we try to find the biggest subset that can be statically chosen and do dynamic index reduction for the remaining equations.

This case is different. The pump model has no dynamic states, so there is no state variable change to be made here, because there are no states at all. Again, I'm not sure how this maps into your algorithms, but as I see it, that's where the difference is.

From what i understand it could be the following: If an MSSS is found that only contains algebraic variables and no natural states, all variables in that set that have been forced to be states will be made algebraic variables, because there is no need in forcing them to be states. Oftentimes the structural singularity of that set will disappear afterwards.

By 'natural states' i mean those which originally appear derived and not those which are forced in by StateSelect.always or StateSelect.prefer.

[casella]

Replying to Karim.Abdelhak:

Yes this gets quite convoluted. The reason we want to differentiate w.r.t. the candidates is that we look for solvability.

OK. Computing the Jacobian of the constraint equations and checking that it is regular is one way to do it. However, if the new states could be computed explicitly from the old states via a sequence of explicit assignments (-> strictly lower triangular BLT), wouldn't that be enough to declare the system solvable, even though the assignments involve function for which we don't have analytic partial derivatives? I guess that would be the case here.

From what i understand it could be the following: If an MSSS is found that only contains algebraic variables and no natural states, all variables in that set that have been forced to be states will be made algebraic variables, because there is no need in forcing them to be states. Oftentimes the structural singularity of that set will disappear afterwards.

By 'natural states' i mean those which originally appear derived and not those which are forced in by StateSelect.always or StateSelect.prefer.

Sounds good. Does this solve our problem?

[casella]

Based on e-mail discussion with Karim, it seems that there are problems in the backend due to the differentiation of expressions containing the homotopy() operator.

The simplified model in homotopy(actual, simplified) is meant to help solving the initialization problem, particularly when it involves large systems of nonlinear equations. It plays absolutely no role during dynamic simulation.

Hence, when dealing with index reduction and state selection, all instances of homotopy(actual,simplified) should be replaced by the actual argument. Ditto for the BLT analysis of the dynamic equations, both in ODE mode and DAE mode.

The only place where homotopy() is relevant is the solution of the initialization problem.

[vitalij]

Replying to casella:

do I understand it right? Based on the solution from #5170 homotopy should be inlined with actual inside the special case inside Index reduxtion, right?

[Karim.Abdelhak]

Replying to vitalij:

Replying to casella:

do I understand it right? Based on the solution from #5170 homotopy should be inlined with actual inside the special case inside Index reduxtion, right?

I tried updating the general handling of the homotopy operator regarding differentiation and how the adjacency/incidence matrix entries are generated. I basically tried to only look at the actual expression and not the simplified one. Furthermore i tried removing the homoptopy operator when differentiating, such that only the differentiated version of the actual expression remains:

homotopy(actual, simplified) --d/dt--> der(actual)

is this the correct approach or should it be like

homotopy(actual, simplified) --d/dt--> homotopy(der(actual), der(simplified))

This should be generally faster than inlineing while producing the same results, unfortunately it does not seem to have an impact on the HeatingSystem we are looking at.

The other system @Francesco send me FlexiCaL.FeedWater.Tests.Tests_FWPHs_Accurate.Tests_with_turbine.Test_FWPHs_14_13_12_11 seems to run with the first approach. Since it is crazy verbose i won't try to analyze any of the results ... i will try to push it and leave that to Francesco.

Regarding the HeatingSystem, this is what my current implementation reports:

### The Equation ### 
0.0 = if tank.regularFlow[2] then tank.ports[2].p - homotopy(tank.vessel_ps_static[2] + 0.5 * tank.portAreas[2] ^ (-2.0) * Modelica.Fluid.Utilities.regSquare2(-m_flow, tank.m_flow_small, (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2], (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.ports_penetration[2] * tank.heatTransfer.states[1].d), false, 1.0), tank.vessel_ps_static[2]) else if tank.inFlow[2] then tank.ports[2].p - tank.vessel_ps_static[2] else -m_flow

--- could not be differentiated for artificial variable ---
 pump.medium.h.

### The Equation ### 
0.0 = if tank.regularFlow[2] then tank.ports[2].p - homotopy(tank.vessel_ps_static[2] + 0.5 * tank.portAreas[2] ^ (-2.0) * Modelica.Fluid.Utilities.regSquare2(-m_flow, tank.m_flow_small, (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2], (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.ports_penetration[2] * tank.heatTransfer.states[1].d), false, 1.0), tank.vessel_ps_static[2]) else if tank.inFlow[2] then tank.ports[2].p - tank.vessel_ps_static[2] else -m_flow

--- could not be differentiated for artificial variable ---
 pump.medium.p.

### The Equation ### 
heater.flowModel.mus[1] = Modelica.Media.Water.IF97_Utilities.dynamicViscosity(heater.statesFM[1].d, heater.statesFM[1].T, pump.medium.p, 0)

--- could not be differentiated for artificial variable ---
 pump.medium.p.

### The Equation ### 
0.0 = if tank.regularFlow[2] then tank.ports[2].p - homotopy(tank.vessel_ps_static[2] + 0.5 * tank.portAreas[2] ^ (-2.0) * Modelica.Fluid.Utilities.regSquare2(-m_flow, tank.m_flow_small, (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2], (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.ports_penetration[2] * tank.heatTransfer.states[1].d), false, 1.0), tank.vessel_ps_static[2]) else if tank.inFlow[2] then tank.ports[2].p - tank.vessel_ps_static[2] else -m_flow

--- could not be differentiated for artificial variable ---
 tank.level.

Equation 1,2 and 4 can't be differentiated, because the homotopy operator contains another function call Modelica.Fluid.Utilities.regSquare2, of which the derivative seems to be unknown. Equation 3 can't be differentiated because the function call Modelica.Media.Water.IF97_Utilities.dynamicViscosity itself seems to have no derivative defined.

This results in revoking the state status of pump.medium.p, pump.medium.h, tank.level which furthermore results in a static state selection. This is actually close to what dymola does, dymola only revokes pump.medium.p and pump.medium.h.
This breaks some of the other tests in the testsuite, but maybe it is still the way to go. Does this seem reasonable?

[Karim.Abdelhak]

The first update regarding homotopy is pushed. As already mentioned it does not fix this particular issue, but another model @Francesco sent me:

FlexiCaL.FeedWater.Tests.Tests_FWPHs_Accurate.Tests_with_turbine.Test_FWPHs_14_13_12_11

Please tell me how the tests went and if the results can be verified!

[casella]

Replying to Karim.Abdelhak:

I tried updating the general handling of the homotopy operator regarding differentiation and how the adjacency/incidence matrix entries are generated. I basically tried to only look at the actual expression and not the simplified one.

That is one option, if you want to avoid having multiple variants of the same equation in memory to save space.

Furthermore i tried removing the homoptopy operator when differentiating, such that only the differentiated version of the actual expression remains:

homotopy(actual, simplified) --d/dt--> der(actual)

is this the correct approach

Of course, it is consistent with only looking at the actual expression when generating incidence/adjacency information. Basically, the simplified argument is only relevant in these two cases:

• when generating the overall set of initialization equations for lambda = 0 ; in this case, actual should be inlined, at least conceptually
• when generating the overall set of initialization equation using homotopy; in this case, homotopy should be inlined as lambda*actual + (1-lambda)*simplified

homotopy(actual, simplified) --d/dt--> homotopy(der(actual), der(simplified))

If you are differentiating equations for the sake of state variable selection and index reduction, der(simplified) is irrelevant. In principle, Pantelides' algorithm can also be used to determine how many initial equations you must add to get a closed initialization probelm, so one could also check what happens with the simplifie constraint equations. However, I would say that the simplifications in the equations for homotopy should not change the structure of the system in terms of state selection, so I would suggest we assume a priori that the state selection performed with actual still leads to a solvable system when lambda = 0. If not, we'll get a runtime error, so be it.

This should be generally faster than inlining while producing the same results,

I guess so.

Unfortunately it does not seem to have an impact on the HeatingSystem we are looking at.

Too bad.

The other system @Francesco send me FlexiCaL.FeedWater.Tests.Tests_FWPHs_Accurate.Tests_with_turbine.Test_FWPHs_14_13_12_11 seems to run with the first approach.

That's very good news, I'm eager to try that out as soon as your changes have been pushed to the latest nighly.

Regarding the HeatingSystem, this is what my current implementation reports:

### The Equation ### 
0.0 = if tank.regularFlow[2] then tank.ports[2].p - homotopy(tank.vessel_ps_static[2] + 0.5 * tank.portAreas[2] ^ (-2.0) * Modelica.Fluid.Utilities.regSquare2(-m_flow, tank.m_flow_small, (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2], (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.ports_penetration[2] * tank.heatTransfer.states[1].d), false, 1.0), tank.vessel_ps_static[2]) else if tank.inFlow[2] then tank.ports[2].p - tank.vessel_ps_static[2] else -m_flow

--- could not be differentiated for artificial variable ---
 pump.medium.h.

### The Equation ### 
0.0 = if tank.regularFlow[2] then tank.ports[2].p - homotopy(tank.vessel_ps_static[2] + 0.5 * tank.portAreas[2] ^ (-2.0) * Modelica.Fluid.Utilities.regSquare2(-m_flow, tank.m_flow_small, (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2], (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.ports_penetration[2] * tank.heatTransfer.states[1].d), false, 1.0), tank.vessel_ps_static[2]) else if tank.inFlow[2] then tank.ports[2].p - tank.vessel_ps_static[2] else -m_flow

--- could not be differentiated for artificial variable ---
 pump.medium.p.

### The Equation ### 
heater.flowModel.mus[1] = Modelica.Media.Water.IF97_Utilities.dynamicViscosity(heater.statesFM[1].d, heater.statesFM[1].T, pump.medium.p, 0)

--- could not be differentiated for artificial variable ---
 pump.medium.p.

### The Equation ### 
0.0 = if tank.regularFlow[2] then tank.ports[2].p - homotopy(tank.vessel_ps_static[2] + 0.5 * tank.portAreas[2] ^ (-2.0) * Modelica.Fluid.Utilities.regSquare2(-m_flow, tank.m_flow_small, (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2], (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.ports_penetration[2] * tank.heatTransfer.states[1].d), false, 1.0), tank.vessel_ps_static[2]) else if tank.inFlow[2] then tank.ports[2].p - tank.vessel_ps_static[2] else -m_flow

--- could not be differentiated for artificial variable ---
 tank.level.

Equation 1,2 and 4 can't be differentiated, because the homotopy operator contains another function call Modelica.Fluid.Utilities.regSquare2, of which the derivative seems to be unknown. Equation 3 can't be differentiated because the function call Modelica.Media.Water.IF97_Utilities.dynamicViscosity itself seems to have no derivative defined.

This results in revoking the state status of pump.medium.p, pump.medium.h, tank.level which furthermore results in a static state selection. This is actually close to what dymola does, dymola only revokes pump.medium.p and pump.medium.h.

I'm fine with revoking the status of pump.medium.p and pump.medium.h, since the pump has no dynamic behaviour. Revoking tank.level instead is not a good idea, that is the best choice of state we have.

(One of) the reason(s) tank.level is revoked potential state status is the lack of a derivative for regRoot2. That function is a one-liner, which should be easily inlined and differentiated, but unfortunately it calls a utility function which is not. Anyway, Dymola can actually differentiate it, and I guess we should, too, because regRoot2 is widely used in Modelica.Fluid.

For example, if I run this test model

model M
  Real y,q;
equation 
  der(y) = q;
  y = Modelica.Fluid.Utilities.regRoot2(time);
end M;

in Dymola, I get

Differentiated the equation
y = smooth(2, smooth(2, (if noEvent(time >= 0.01) then sqrt(time) else (if noEvent(time <= -0.01) then -sqrt(abs(time)) else Modelica.Fluid.Utilities.regRoot2.regRoot2_utility (time, 0.01, 1, 1, false, 1)))));

giving
der(y) = smooth(1, smooth(1, (if noEvent(time >= 0.01) then 0.5/sqrt(time) else (if noEvent(time <= -0.01) then -0.5*noEvent((if time > 0 then 1 else -1))/sqrt(abs(time)) else Modelica.Fluid.Utilities.regRoot2.regRoot2_utility:der (time, 0.01, 1, 1, false, 1, 1.0, 0.0, 0, 0, 0)))));

[Karim.Abdelhak]

Replying to casella:

For example, if I run this test model

model M
  Real y,q;
equation 
  der(y) = q;
  y = Modelica.Fluid.Utilities.regRoot2(time);
end M;

in Dymola, I get

Differentiated the equation
y = smooth(2, smooth(2, (if noEvent(time >= 0.01) then sqrt(time) else (if noEvent(time <= -0.01) then -sqrt(abs(time)) else Modelica.Fluid.Utilities.regRoot2.regRoot2_utility (time, 0.01, 1, 1, false, 1)))));

giving
der(y) = smooth(1, smooth(1, (if noEvent(time >= 0.01) then 0.5/sqrt(time) else (if noEvent(time <= -0.01) then -0.5*noEvent((if time > 0 then 1 else -1))/sqrt(abs(time)) else Modelica.Fluid.Utilities.regRoot2.regRoot2_utility:der (time, 0.01, 1, 1, false, 1, 1.0, 0.0, 0, 0, 0)))));

We actually get something similar in OpenModelica using -d=debugDifferentiation (stripping some unnecessary verbose output):

### differentiateEquationTime
y = Modelica.Fluid.Utilities.regRoot2(time, 0.01, 1.0, 1.0, false, 1.0) w.r.t. time

### Result of differentiateEquationTime
 --> der(y) = smooth(1, if noEvent(time >= 0.01) then 0.5 / sqrt(time) else if noEvent(time <= -0.01) then (-0.5) * sign(time) / sqrt(abs(time)) else $DER$$PModelica$PFluid$PUtilities$PregRoot2$PregRoot2_utility(time, 0.01, 1.0, 1.0, false, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0))

I am not quite sure if the results are equal. We seem to strip one call of smooth, i guess that is done by some expression simplification module. The last part is just that dymola seems to split up the sign operator and we don't, and the expression for the function derivative is different:

Modelica.Fluid.Utilities.regRoot2.regRoot2_utility:der
$DER$$PModelica$PFluid$PUtilities$PregRoot2$PregRoot2_utility

The only thing i could think of is the missing smooth operator, but i guess stating twice that it is once differentiable should not make any difference and be the same as stating it once.

I think i have to find out why the differentiation seems to fail for the HeatingSystem and not in this minimal case, maybe it is because of the nested function calls.

[casella]

Replying to Karim.Abdelhak:

We actually get something similar in OpenModelica using -d=debugDifferentiation (stripping some unnecessary verbose output):

### differentiateEquationTime
y = Modelica.Fluid.Utilities.regRoot2(time, 0.01, 1.0, 1.0, false, 1.0) w.r.t. time

### Result of differentiateEquationTime
 --> der(y) = smooth(1, if noEvent(time >= 0.01) then 0.5 / sqrt(time) else if noEvent(time <= -0.01) then (-0.5) * sign(time) / sqrt(abs(time)) else $DER$$PModelica$PFluid$PUtilities$PregRoot2$PregRoot2_utility(time, 0.01, 1.0, 1.0, false, 1.0, 1.0, 0.0, 0.0, 0.0, 0.0))

Aha, nice. I should have tried this out myself :)

I am not quite sure if the results are equal. We seem to strip one call of smooth

That's no big deal

, i guess that is done by some expression simplification module. The last part is just that dymola seems to split up the sign operator and we don't, and the expression for the function derivative is different:

Modelica.Fluid.Utilities.regRoot2.regRoot2_utility:der
$DER$$PModelica$PFluid$PUtilities$PregRoot2$PregRoot2_utility

I guess this refers to some automatically generated derivative function, since regRoot2_utility has no derivative annotation. Of course the name of such a function is completely arbitrary.

The only thing i could think of is the missing smooth operator, but i guess stating twice that it is once differentiable should not make any difference and be the same as stating it once.

Absolutely.

I think i have to find out why the differentiation seems to fail for the HeatingSystem and not in this minimal case, maybe it is because of the nested function calls.

Yeah, I also thought about the nested function calls, please do check that.

[Karim.Abdelhak]

I found something quite disturbing... please consider following model:

model M
  Real y(stateSelect=StateSelect.prefer);
  Real x, x2;
equation
  x2 = sin(time);
  der(x) = y^2;
  y = Modelica.Fluid.Utilities.regRoot2(x, x2);
end M;

This runs perfectly fine in OpenModelica but changing the third equation to a residual by hand breaks it:

model M
  Real y(stateSelect=StateSelect.prefer);
  Real x, x2;
equation
  x2 = sin(time);
  der(x) = y^2;
  0 = Modelica.Fluid.Utilities.regRoot2(x, x2) - y;
end M;

Also changing the first equation from x2 = sin(time) to x2 = cos(time) somehow avoids the error (???).

I had a deeper look into this and realized that only differentiation w.r.t. to a specific variable does not work, differentiating w.r.t. time is no problem. The rhs and lhs are processed individually, therefore allowing a matching for y in the first case, whereas in the latter case the differentiation of the function and therefore for the whole rhs fails. The result is that no candidate can be solved in this equation at all.

I think i have to dive deep into the differentiation module and remove that error. If it is differentiable by time it is also differentiable by any variable.

[casella]

Replying to Karim.Abdelhak:

If it is differentiable by time it is also differentiable by any variable.

Unless the differentiation involves functions with derivative annotation, which only provide the total time derivative, not the individual ones.

Of course one may call the derivative function N times, each time setting one input derivative to 1 and the others to 0. This actually gives you the partial derivatives. But that's not very efficient from a computational standpoint.

[Karim.Abdelhak]

Replying to casella:

Of course one may call the derivative function N times, each time setting one input derivative to 1 and the others to 0. This actually gives you the partial derivatives. But that's not very efficient from a computational standpoint.

For the generation of symbolic jacobians we oftentimes have no choice but to do it. We kind of optimized it there so i think we are safe to say "if we can derive by time we can derive by a variable".

In our case it is just a general solvability check where the actual result does not matter - only if it is linear, nonlinear, constant or zero. And a derivative of zero is a great difference to the solvability check failing completely. But i fixed that - above model should now work in any constellation.

For our HeatingSystem the change results in only reverting pump.medium.p due to the dynamicViscosity function. This unfortunately does not suffice and the model still fails, but i found following equation:

pump.Hb_flow = m_flow * (smooth(0, tank.medium.h) - smooth(0, pump.medium.h))

it clearly states that pump.medium.h should not be differentiated, so i implemented a check for such smooth(0, ...) calls on hard artificial states and revert them. Works fine for this model, let's see what the rest of the testsuite looks like.

For now it only checks if the artifical state should not be differentiated on it's own like smooth(0, pump.medium.h), should this also trigger if it is contained in any expression like
smooth(0, k*pump.medium.h^2) or any other expression you could think of containing the artificial state?

Any thoughts on this?

[casella]

Replying to Karim.Abdelhak:

This unfortunately does not suffice and the model still fails

Boy, this is clearly going to be the longest ticket ever. Are we aiming at comment:100? :)

but i found following equation:

pump.Hb_flow = m_flow * (smooth(0, tank.medium.h) - smooth(0, pump.medium.h))

it clearly states that pump.medium.h should not be differentiated, so i implemented a check for such smooth(0, ...) calls on hard artificial states and revert them. Works fine for this model, let's see what the rest of the testsuite looks like.

Hmm, this is quite weird. Were does this equation come from? Is it from the source code, or was it derived by the backend?

Besides, this is the definition of smooth in the Specification:

If p>=0 smooth(p,expr) returns expr and states that expr is p times continuously differentiable, i.e.: expr is continuous in all real variables appearing in the expression and all partial derivatives with respect to all appearing real variables exist and are continuous up to order p

My interpretation is that the expression tank.medium.h is the same as x, i.e., just one single variable. Stating smooth(0,x) seems quite a stupid thing to do, since x is obviously an analytical function and can be continuously differentiated as many time as you want. The first derivative is one, and all higher-order derivatives are zero.

Stated in other words, the smooth(<expr>) operator states that <expr> is smooth with respect to its variables, not with respect to time.

Do I miss something?

For now it only checks if the artifical state should not be differentiated on it's own like smooth(0, pump.medium.h), should this also trigger if it is contained in any expression like
smooth(0, k*pump.medium.h^2) or any other expression you could think of containing the artificial state?

See above. I don't see the point of applying smooth to polynomials, that are obviously infinitely times differentiable in closed form.

[Karim.Abdelhak]

Replying to casella:

Replying to Karim.Abdelhak:

This unfortunately does not suffice and the model still fails

Boy, this is clearly going to be the longest ticket ever. Are we aiming at comment:100? :)

Looks like it, but with the current fix (smooth operator testing) everything works.

Hmm, this is quite weird. Were does this equation come from? Is it from the source code, or was it derived by the backend?

The original equation seems to be

 pump.Hb_flow = pump.port_a.m_flow * smooth(0, tank.ports[2].h_outflow) + pump.port_b.m_flow * smooth(0, pump.port_b.h_outflow)

and the alias module seems to replace some expressions so that we end up with the other equation.

Besides, this is the definition of smooth in the Specification:

If p>=0 smooth(p,expr) returns expr and states that expr is p times continuously differentiable, i.e.: expr is continuous in all real variables appearing in the expression and all partial derivatives with respect to all appearing real variables exist and are continuous up to order p

I mean i just went through kind of the whole system and tried to determine why dymola decides to revoke this variables state status and this seems to be the only equation where one could argue. You are right that it does not make sense on polynomials, but this should not end up in index reduction.

As i said, everything seems to work with my fix - should i consult some other developers to discuss about this?

​PR419

[casella]

Replying to Karim.Abdelhak:

The original equation seems to be

pump.Hb_flow = pump.port_a.m_flow * smooth(0, tank.ports[2].h_outflow) +
               pump.port_b.m_flow * smooth(0, pump.port_b.h_outflow)

Aha, wait a minute. The original equations are found in ​PartialPump:

    Hb_flow = port_a.m_flow*actualStream(port_a.h_outflow) +
              port_b.m_flow*actualStream(port_b.h_outflow);

According to ​Section 15.3 of the Specification:

actualStream(c.h_outflow) = if c.m_flow > 0 
  then inStream(c.h_outflow)
  else c.h_outflow;

I guess the backend somehow infers that the mass flow rate has a fixed direction, so the conditional equation is gone. However, I have no idea why smooth should be added at all. It's good that everything works with your PR, but I am not sure that relying on a dubious interpetation of a smooth operator added by the backend for some unexplained reason is a good idea.

In fact, there could be an explanation for the conditional equation being simplified away, because the pump model has allowFlowReversal = false, but I still don't understand why smooth should be added at all.

Do I miss something?

[Karim.Abdelhak]

Replying to casella:

I guess the backend somehow infers that the mass flow rate has a fixed direction, so the conditional equation is gone. However, I have no idea why smooth should be added at all. It's good that everything works with your PR, but I am not sure that relying on a dubious interpetation of a smooth operator added by the backend for some unexplained reason is a good idea.

Right below the specification of actualStream operator we can find

In the case of equation (1), although the actualStream() operator is discontinuous, the product with the flow variable is not, because actualStream() is discontinuous when the flow is zero by construction. Therefore, a tool might infer that the expression is smooth(0, ...) automatically, and decide whether or not to generate an event. [...]

That probably is the reason for the smooth operator. I will consult some others about this.

Following equations are the ones on which i can decide to revoke the state status of pump.medium.h, maybe you got an idea why dymola decides to do it?

1/1 (1): Modelica.Media.Water.IF97_Utilities.rho_props_ph(pump.medium.p, pump.medium.h, Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(pump.medium.p, pump.medium.h, pump.heatTransfer.states[1].phase, 0)) - pump.rho=0   [dynamic |0|0|0|0|]
2/2 (2): m_flow * (smooth(0, tank.medium.h) - smooth(0, pump.medium.h)) - pump.Hb_flow=0   [dynamic |0|0|0|0|]
3/3 (3): Modelica.Media.Water.IF97_Utilities.rho_props_ph(pump.medium.p, pump.medium.h, Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(pump.medium.p, pump.medium.h, 0, 0)) - heater.statesFM[1].d= 0   [dynamic |0|0|0|0|]
4/4 (4): Modelica.Media.Water.IF97_Utilities.T_props_ph(pump.medium.p, pump.medium.h, Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(pump.medium.p, pump.medium.h, 0, 0)) - heater.statesFM[1].T= 0   [dynamic |0|0|0|0|    ]

They are transformed to residual form.

[casella]

Replying to Karim.Abdelhak:

Right below the specification of actualStream operator we can find

In the case of equation (1), although the actualStream() operator is discontinuous, the product with the flow variable is not, because actualStream() is discontinuous when the flow is zero by construction. Therefore, a tool might infer that the expression is smooth(0, ...) automatically, and decide whether or not to generate an event. [...]

Yeah, sorry, I took part in writing that myself, I should have thought about it. The smooth operator should then be applied to the entire expression, i.e.

pump.Hb_flow = smooth(0, pump.port_a.m_flow * tank.ports[2].h_outflow) +
               smooth(0, pump.port_b.m_flow * pump.port_b.h_outflow)

however, this really only makes sense if actualStream is replaced by a conditional expression; if that expression is simplified because of assumptions on the mass flow rate sign (which is good!), then you get a plain product, which is infinitely times differentiable, so there is no longer any reason to apply smooth to it. Also, I don't think it makes any sense to symbolically pull out pump.port_a.m_flow from the smooth expression, but non the enthalpy term.

Following equations are the ones on which i can decide to revoke the state status of pump.medium.h, maybe you got an idea why dymola decides to do it?

1/1 (1): Modelica.Media.Water.IF97_Utilities.rho_props_ph(pump.medium.p, pump.medium.h, Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(pump.medium.p, pump.medium.h, pump.heatTransfer.states[1].phase, 0)) - pump.rho=0   [dynamic |0|0|0|0|]
2/2 (2): m_flow * (smooth(0, tank.medium.h) - smooth(0, pump.medium.h)) - pump.Hb_flow=0   [dynamic |0|0|0|0|]
3/3 (3): Modelica.Media.Water.IF97_Utilities.rho_props_ph(pump.medium.p, pump.medium.h, Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(pump.medium.p, pump.medium.h, 0, 0)) - heater.statesFM[1].d= 0   [dynamic |0|0|0|0|]
4/4 (4): Modelica.Media.Water.IF97_Utilities.T_props_ph(pump.medium.p, pump.medium.h, Modelica.Media.Water.IF97_Utilities.waterBaseProp_ph(pump.medium.p, pump.medium.h, 0, 0)) - heater.statesFM[1].T= 0   [dynamic |0|0|0|0|    ]

I have no idea, and in fact I'm not sure that Dymola takes the decision based on exactly these equations, which should in principle be all differentiable. It could even be that Dymola does not infer that the mass flow rate inside actualStream is positive, so it will keep the conditional expression inside smooth, and that would be enough to revoke the potential state status.

According to my understanding, the reason why tank.medium.h should be preferred over pump.medium.h in this case (whereby the conditional equation has been removed) is that solving the equations w.r.t. pump.medium.h and computing all the required derivatives will be in general more involved than solving them for tank.medium.h. This is because tank is the component with dynamic balance equations, so using variables from another component as states will involve a lot more algebraic equations and variables in the process.

In any case, the choice of tank.medium.h can lead to a static state selection, while apparently the choice of pump.medium.h leads to a dynamic state selection, which was the original reason for opening this ticket. Isn't this enough cause to prefer tank.medium.h?

[Karim.Abdelhak]

Replying to casella:

In any case, the choice of tank.medium.h can lead to a static state selection, while apparently the choice of pump.medium.h leads to a dynamic state selection, which was the original reason for opening this ticket. Isn't this enough cause to prefer tank.medium.h?

In theory it is, but the problem is, that keeping pump.medium.h as potential state forces the equation

1/1 (1): 0.0 = if tank.regularFlow[2] then tank.ports[2].p - homotopy(tank.vessel_ps_static[2] + 0.5 * tank.portAreas[2] ^ (-2.0) * smooth(2, if noEvent((-m_flow) >= tank.m_flow_small) then (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2] * m_flow ^ 2.0 else if noEvent((-m_flow) <= (-tank.m_flow_small)) then m_flow ^ 2.0 * (-1.0 + (tank.portAreas[2] / tank.vesselArea) ^ 2.0 - tank.portsData[2].zeta_out) / (tank.heatTransfer.states[1].d * tank.ports_penetration[2]) else if noEvent((-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2] >= (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[2])) then Modelica.Fluid.Utilities.regSquare2.regSquare2_utility(-m_flow, tank.m_flow_small, (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2], (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[2]), false, 1.0) else Modelica.Fluid.Utilities.regSquare2.regSquare2_utility(m_flow, tank.m_flow_small, (1.0 + tank.portsData[2].zeta_out - (tank.portAreas[2] / tank.vesselArea) ^ 2.0) / (tank.heatTransfer.states[1].d * tank.ports_penetration[2]), (-1.0 + tank.portsData[2].zeta_in + (tank.portAreas[2] / tank.vesselArea) ^ 2.0) * tank.ports_penetration[2] / tank.portInDensities[2], false, 1.0)), tank.vessel_ps_static[2]) else if tank.inFlow[2] then tank.ports[2].p - tank.vessel_ps_static[2] else -m_flow [unknown |0|0|0|0|]

into index reduction. And this one causes all the problems since it can't really be solved for anything explicitly.

To be exact the following MSSS shrunk from

##############--MSSS--##############                                                                               
Indices of constraint equations: 32 3 33 5 39 43 50 51 80 54 49 46 53 29 23 25 11 55 99 104 77 78 88 98 74 86 73 10 5 76  

to

##############--MSSS--##############
Indices of constraint equations: 32 3 33 5

because we removed pump.medium.h before even starting index reduction. As a simple example consider the following (abstract) model where someWierdFunction is supposed to represent a function that causes problems for index reduction:

model M
  Real x, y, z(stateSelect=StateSelect.prefer);
equation
  der(y) = sin(time);
  someWierdFunction(x) = 0.0;
  x = z + y;
end M;

If we analyze this with pantelides algorithm we come to the conclusion, that we have to differentiate the second and third equation, because those two have only one unknown x. This results in someWierdFunction being part of an MSSS. If we revoked z beforehand, this would not have happened and the causalization would be quite straight foreward in x being implicitly solved in equation (2) and z being solved in equation (3).

Somehow dymola seems to manage to find this conclusion.

[Karim.Abdelhak]

I thought once again about the smooth operator. I will address this in today's meeting since we seem to disagree on this, but i want to shortly summarize what i think about this.

As you quoted from the specification, any expression x inside smooth(0, x) is zero times continuously differentiable w.r.t. to all real expressions contained in it. In our case we check for partial derivatives w.r.t. all state candidates to check for solvability. Furthermore we would need those partial derivatives for the symbolic jacobian. What i implemented is a "sanity check" for all artificial variables (forced in by stateSelect and not appearing derived before index reduction) and see if they would cause problems for index reduction. Since smooth(0, pump.medium.h) implies that it is zero times differentiable by all real expressions, that includes pump.medium.h itself and therefore it can't be differentiated w.r.t. to it -> the sanity check fails and revokes pump.medium.h. I think this is the right approach since it should never end up in index reduction because we cannot derive it partially w.r.t. to all state candidates.

[casella]

Summary of today's webmeeting

• in this specific case smooth should be gone because there are no longer any conditional equations due to min/max constraint on the connector flow variables
• there are probably other issues in the way the frontend handles inStream, we will open a ticket about that for @perost to take care of
• once the frontend issues are solved, PR 419 should be ok, i.e., if smooth(0, if c.m_flow > 0 then something else something else actually shows up in the flattened code (as would happen in HeatingSystem if the pump had allowFlowReversal = true), the expression will indeed be not differentiable and should cause the revoking of the variables in something and something else from the list of potential states
• in general the front-end should avoid generating expressions which are obviously infinitely times differentiable inside smooth(0, ...)
• in general, the backend can indeed assume that smooth(0,<any_expression>) is not continuously differentiable, without looking into the expression itself
• if an end user writes smooth(0, <expr>) and expr is infinitely times diffentiable, too bad for him/her.
• it goes without saying that OMC should not generate expressions with smooth that are obviously infinitely times differentiable, as it currently happens (see, e.g., comment:69)

[vitalij]

Replying to casella:

• in general, the backend can indeed assume that smooth(0,<any_expression>) is not continuously differentiable, without looking into the expression itself
• if an end user writes smooth(0, <expr>) and expr is infinitely times diffentiable, too bad for him/her.
• it goes without saying that OMC should not generate expressions with smooth that are obviously infinitely times differentiable, as it currently happens (see, e.g., comment:69)

Thanks for the summary!
I did not quite understand the topic with the smooth. What is the impact of other piecewise function like abs? Also function which are not continuously differentiable but piecewise differentiable. I thought the idea of ​​smooth was for the runtime and e.g. for the stepsize or/and order control or how was that again?

[casella]

Replying to vitalij:

I did not quite understand the topic with the smooth. What is the impact of other piecewise function like abs?
Also function which are not continuously differentiable but piecewise differentiable.

Of course they are also non continuously diffentiable, so this has an impact on Pantelides' algorithm.

I thought the idea of ​​smooth was for the runtime and e.g. for the stepsize or/and order control or how was that again?

That is one possible use of the operator, but of course the information is relevant also when differentiating for index reduction. This also holds for the smoothOrder function annotation.

[casella]

See #5642

[vitalij]

For me it's not trivial to understand, because of the Picard–Lindelöf theorem say we need
Lipschitz continuity and not continuously differentiable. A piecewise differentiable function can sensing Lipschitz continuity.

[casella]

Replying to vitalij:

For me it's not trivial to understand, because of the Picard–Lindelöf theorem say we need
Lipschitz continuity and not continuously differentiable.

This regards the existence of the solution of the differential equations.

The discussion here is about index reduction and state variable change, which has different requirements. I'm not 100% sure continuous differentiability is necessary, but I understand it's safe to require it for Pantelides' algorithm. @karim, can you comment on that?

[casella]

Unfortunately ​PR #419 caused some ​regressions, mostly concerning mechanical systems. These are not related to the actualStream fix #5642, so they should be investigated further.

A representative case from MSL is ​Modelica.Mechanics.Translational.Examples.HeatLosses. The backend reports:

Warning: Some equations could not be differentiated for following variables having attribute stateSelect=StateSelect.prefer. They will be treated as if they had stateSelect=StateSelect.default
========================================
1: damper.v_rel
2: elastoGap.v_rel
3: springDamper.v_rel
4: springDamper1.v_rel

which later causes a failure of Pantelides' algorithm.

This is the static state selection performed by Dymola:

damper.s_rel
damper.v_rel
massWithStopAndFriction.s
massWithStopAndFriction.v
springDamper.s_rel
springDamper.v_rel
springDamper1.s_rel
springDamper1.v_rel

I guess the currently implemented criterion for excluding candidate states is too strict for some reason.

[casella]

Notice that in the ​Modelica.Mechanics.Translational.Examples.HeatLosses case, the simulation is carried out successfully, but the verification, which was previously successful, now fails.

If you check the ​verification report, you will see the same kind of numerical errors that plague #5328. I guess this is by no means a coincidence :)

[casella]

I checked the -d=bltdump log. I found messages such as:

### The Equation ###  
springDamper1.v_rel = der(springDamper1.flange_b.s) - brake.v 
 
--- could not be differentiated for artificial variable --- 
 springDamper1.v_rel. 

Why not? The partial derivative springDamper1.v_rel w.r.t. springDamper1.v_rel is one. Is this because it can't handle der(der(springDamper1.flange_b.s))?

[Karim.Abdelhak]

Replying to casella:

I checked the -d=bltdump log. I found messages such as:

### The Equation ###  
springDamper1.v_rel = der(springDamper1.flange_b.s) - brake.v 
 
--- could not be differentiated for artificial variable --- 
 springDamper1.v_rel. 

Why not? The partial derivative springDamper1.v_rel w.r.t. springDamper1.v_rel is one. Is this because it can't handle der(der(springDamper1.flange_b.s))?

I think the differentiation and adjacency matrix module is a little bit buggy and fails on more cases than it should. In this case the der call is not handled properly. I will have a look into that!

[casella]

Excellent! I hope this solves some issues in the testsuite :)

[Karim.Abdelhak]

Quite possibly fixed with ​PR457.
I added replacement of the derivative operator because the differentiation module apparently can't handle it.

[casella]

Just started the ​newInst Jenkins job, let's see what happened when it is completed.

[Karim.Abdelhak]

Replying to casella:

Just started the ​newInst Jenkins job, let's see what happened when it is completed.

The ​results look good so far ... i guess? Is there something left to be fixed?

And are the two regressions crucial?

[casella]

Those two regressions are bogus, they come and go each time the testsuite is run, no idea why.

My only concern is that ​the previous commit had 66 regressions, while ​your last fix only shows 55 improvements.

For example, most variants of ModelicaTest.Fluid.TestExamplesVariants.BranchingDynamicPipes_MomentumSteadyState got broken first, but they were not restored later. I hope this is due to the smooth issue #5642, which is still affecting Modelica.Fluid.Examples.HeatingSystem. Let's wait until that is fixed before closing this ticket.

Unfortunately, #5328 is still valid. You may want to have another look at that.

[Karim.Abdelhak]

Replying to casella:

For example, most variants of ModelicaTest.Fluid.TestExamplesVariants.BranchingDynamicPipes_MomentumSteadyState got broken first, but they were not restored later. I hope this is due to the smooth issue #5642, which is still affecting Modelica.Fluid.Examples.HeatingSystem. Let's wait until that is fixed before closing this ticket.

I will check that, i want to implement the full sanity check anyways, even if the frontend part of this is not fixed yet.

Unfortunately, #5328 is still valid. You may want to have another look at that.

Yea, maybe it is connected to this, but it really is quite the mistery. Probably also something in the differentiation module such that the jacobian is slightly wrong. Andreas and I are working on the delay operator differentiation, that actually produces wrong jacobians so it could be the origin of the problem.

[casella]

Replying to Karim.Abdelhak:

Unfortunately, #5328 is still valid. You may want to have another look at that.

Yea, maybe it is connected to this, but it really is quite the mistery. Probably also something in the differentiation module such that the jacobian is slightly wrong. Andreas and I are working on the delay operator differentiation, that actually produces wrong jacobians so it could be the origin of the problem.

There should be no delay operators there.

My recommendation is that you implement a debug flag that compares numerical jacobians with symbolic jacobians - this would probably help finding the bug.

[Karim.Abdelhak]

Replying to casella:

My recommendation is that you implement a debug flag that compares numerical jacobians with symbolic jacobians - this would probably help finding the bug.

Sounds like a general good idea, but it actually isn't that easy to do both at the same time and compare them. I will try to figure something out with Andreas.

[casella]

Replying to Karim.Abdelhak:

Replying to casella:

My recommendation is that you implement a debug flag that compares numerical jacobians with symbolic jacobians - this would probably help finding the bug.

Sounds like a general good idea, but it actually isn't that easy to do both at the same time

This may not be necessary. In fact, we could simply run the test case twice, dump the jacobian with some debug flag, and then parse them in some other tool. It's not something that you do all the time. I understand this is already possible now, maybe we should just stick to that.

and compare them.

This may be somewhat tricky, because of the usual scaling issues. Elements of the Jacobian have widely different dimensions, so 1e4 may be a small error for some elements, and 1e-6 a large error for others. Relative errors don't work reliably for elements which are supposed to be zero, and turn out to be 1e-12 in one case and 2e-12 in the other case, of course that's not a 100% error.

This problem would actually be solved automatically if we dump the scaled jacobians, then you can just compute the difference and look for terms which are not <<< 1

[Karim.Abdelhak]

Replying to casella:

Replying to Karim.Abdelhak:

Replying to casella:

My recommendation is that you implement a debug flag that compares numerical jacobians with symbolic jacobians - this would probably help finding the bug.

Sounds like a general good idea, but it actually isn't that easy to do both at the same time

This may not be necessary. In fact, we could simply run the test case twice, dump the jacobian with some debug flag, and then parse them in some other tool. It's not something that you do all the time. I understand this is already possible now, maybe we should just stick to that.

and compare them.

This may be somewhat tricky, because of the usual scaling issues. Elements of the Jacobian have widely different dimensions, so 1e4 may be a small error for some elements, and 1e-6 a large error for others. Relative errors don't work reliably for elements which are supposed to be zero, and turn out to be 1e-12 in one case and 2e-12 in the other case, of course that's not a 100% error.

This problem would actually be solved automatically if we dump the scaled jacobians, then you can just compute the difference and look for terms which are not <<< 1

Well from my understanding we want to simulate everything numerically and additionally dump the symbolical evaluation at the same time (but don't use it!).

The other option would be to kind of interpolate the results and compare from there, but i don't know how we should interpret that information otherwise. In fact Willi implemented something like this for ADOL-C, maybe we could scavenge that.

[casella]

Replying to Karim.Abdelhak:

Well from my understanding we want to simulate everything numerically and additionally dump the symbolical evaluation at the same time (but don't use it!).

Or the other way round, why not?

The other option would be to kind of interpolate the results and compare from there

I'm not sure what you mean here

[Karim.Abdelhak]

Replying to casella:

Replying to Karim.Abdelhak:

Well from my understanding we want to simulate everything numerically and additionally dump the symbolical evaluation at the same time (but don't use it!).

Or the other way round, why not?

Numerically it works and and symbolically it does not, but maybe the other way around would be useful as well.

The other option would be to kind of interpolate the results and compare from there

I'm not sure what you mean here

Run tests with numerical and symbolical jacobians separately and try to compare afterwards. Since the solver step size might not be the same, the values of the jacobians need to be interpolated such that we can compare at the same points in time. But this seems unnecessarily hard, so probably we won't to that.

[Karim.Abdelhak]

As was to be expected this model failed again after #5642 got fixed. With ​PR478 i found a general fix such that the model runs even when not reverting the artificial states and with using dynamic state selection. It also seemed to fix some other minor issues. Nevertheless i will implement the full sanity check to revoke the artificial variables since a static state selection should be preferred if able to do so. After it is merged and we got information about the regressions we hopefully can close this ticket. For the artificial variable check i will make an enhancement ticket.

[casella]

Replying to Karim.Abdelhak:

After it is merged and we got information about the regressions we hopefully can close this ticket.

We need to get to comment:100 before we can do that.

Here is my contribution to the feat :)

[casella]

I checked the results after PR 478. The HeatingSystem model simulates successfully, but it takes a lot of time (3 mins on Jenkins), and if I run it locally, I still see

 Number of states: 10 (
$STATESET1.x[4],$STATESET1.x[3],$STATESET1.x[2],$STATESET1.x[1],
radiator.mediums[1].p,radiator.mediums[1].h,
pipe.mediums[1].p,pipe.mediums[1].h,pipe.mediums[2].p,pipe.mediums[2].h)

@Karim, I'm looking forward to your final sanity check, I hope it resolves the issue for good :)

[casella]

One more comment for the record. While compiling HeatingSystem with the latest nightly build of OMEdit, I noticed several red lines regarding the C code compilation, e.g.

In file included from Modelica.Fluid.Examples.HeatingSystem_model.h:23:0,
                 from Modelica.Fluid.Examples.HeatingSystem_12jac.c:2:
Modelica.Fluid.Examples.HeatingSystem_12jac.c: In function 'Modelica_Fluid_Examples_HeatingSystem_eqFunction_2021':
Modelica.Fluid.Examples.HeatingSystem_12jac.c:276:61: warning: passing argument 2 of 'Modelica_Media_Common_IF97BaseTwoPhase_copy_to_vars_p' from incompatible pointer type [-Wincompatible-pointer-types]
   Modelica_Media_Common_IF97BaseTwoPhase_copy_to_vars(tmp0, &(jacobian->tmpVars[139] /* $cse81.phase JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[138] /* $cse81.region JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[137] /* $cse81.p JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[136] /* $cse81.T JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[135] /* $cse81.h JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[134] /* $cse81.R JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[133] /* $cse81.cp JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[132] /* $cse81.cv JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[131] /* $cse81.rho JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[130] /* $cse81.s JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[129] /* $cse81.pt JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[128] /* $cse81.pd JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[127] /* $cse81.vt JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[126] /* $cse81.vp JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[125] /* $cse81.x JACOBIAN_DIFF_VAR */), &(jacobian->tmpVars[124] /* $cse81.dpT JACOBIAN_DIFF_VAR */));;

I'm not 100% this is relevant to this ticket, but I noticed that when I compile other examples (e.g. DrumBoiler) that have fixed state selection, this doesn't happen, so I see a chance that it actually is.

[casella]

Releasing 1.14.0 which is stable and has many improvements w.r.t. 1.13.2.

This issue, previously marked as blocker for 1.14.0, is rescheduled to 1.15.0

[casella]

We agreed that the sanity check for "questionable" user input, setting StateSelect.prefer where this is really not needed, would be nice to have, but it's not top priority now.

We'll keep the MSL example as it is so we have it as a test case for this feature.

Tentatively rescheduling to 2.0.0

[casella]

After some discussion on ​MSL issue #3236, I decided to fix the HeatingSystem model, and all similar static models that extend PartialLumpedVolume, by making sure that stateSelect.prefer is not set on the pump states, which is just stupid.

This makes omc 1.14.0-dev.beta3 to compute a static state selection :)

As soon as ​PR 3249 is merged into MSL 4.0.0 we will be able to see the impact on the whole testsuite on our Jenkins reports. If all goes fine, we will also backport it to MSL 3.2.3 maintenance and finally close this ticket.

[casella]

The PRs have been merged and static state selection is now achieved as expected