Mixed under/overdetermined initialization problem in Buildings

[Francesco Casella]

Please check ​Buildings.Experimental.DistrictHeatingCooling.Plants.Validation.Plant_Carnot_T_ClosedLoop. The backend fails

Error: The initialization problem of given system is mixed-determined. It is under- as well as overdetermined and the mixed-determination-index is too high. [index > 10]
Please checkout the option "--maxMixedDeterminedIndex" to simulate with a higher threshold or consider changing some initial equations, fixed variables and start values.

We should figure out where does this issue come from.

[Francesco Casella]

Retargeted to 1.18.0 because of 1.17.0 timed release.

[Francesco Casella]

@Karim, it would be nice to have a bit more output besides the mixed-determination index. Would it be possible to know which variables/equations are involved in the under- or over-determined part of the system?

[Karim Adbdelhak]

The information for this specific model would be:

------------ UNBALANCED INITAL SYSTEM ------------
Unresolved overdetermination (unmatched eqns):
	pla.coo.con.vol.T = 273.15 + (pla.coo.port_b1.h_outflow + (-2501014.5) * pla.coo.con.vol.dynBal.medium.Xi[1]) / (1006.0 * (1.0 - pla.coo.con.vol.dynBal.medium.Xi[1]) + 1860.0 * pla.coo.con.vol.dynBal.medium.Xi[1])

Unresolved underdetermination (unmatched vars):
	pum.filter.x[2]:VARIABLE(flow=false stateSelect=StateSelect.always protected = true )  "Filter states" type: Real [2]
--------------------------------------------------

There could be more information extracted, but i would need to dig a lot deeper. Does this information suffice? I would add it to the output of -d=initialization and extend the error message to refer to this debug flag.

FYI: That specific model is moved to Obsolete in the latest buildings master Buildings.Obsolete.DistrictHeatingCooling.Plants.Validation.Plant_Carnot_T_ClosedLoop.

[Francesco Casella]

@Karim reports that this model ended up in an "obsolete" package in Buildings 8.0.0. I have no idea why. @mwetter, could you comment on that?

Should we spend time trying to get this specific model to work, or shold be cut it out from the m7.0 maintenance branch?

[Michael Wetter]

I think the model is a good test for OpenModelica. From a Modelica point of view, I see nothing wrong with the model. The reason why it was moved to Obsolete is because Buildings.Experimental.DHC replaces these systems with a hydronic network that is better suited for 5th generation district heating and cooling. So we don't want users to imitate the hydraulic configuration used in this package.

[Francesco Casella]

@mwetter, thanks for the feedback.

We'll continue investigating what's wrong with it.

[Francesco Casella]

@Karim, as we discussed yesterday, the current ouptut is a bit terse, because we don't know if the "wrong equation" is really the one that remains unmatched at the end of the process. We should be able to see the entire set of candidates for the role of unmatched equations, not just the one that happens to be so at the end, which is likely subject to some arbitrariness.

Ditto for the variables.

[Karim Adbdelhak]

I managed to extract the information about the full set. I tried to present it as well as possible including the fact that we ended up with one unmatched equation and one unmatched variable. It might be the case that there is another problem during initialization though, we do not traverse further to find all. We just stop at the first and report it. Is this what you wanted @Francesco?

Can be accessed via -d=initialization, i will add a note to the error message.

------------ UNBALANCED INITIAL SYSTEM ------------
The initial system is over- as well as underdetermined and it could not be resolved after 10 iterations.

==== OVERDETERMINATION BY 1 EQUATION(S)
==== UNDERDETERMINATION OF 1 VARIABLE(S)

---- involved set eqns (26/25):
  2(2):	algorithm
  pulse.count := integer((time - pulse.startTime) / pulse.period);
  pulse.T_start := pulse.startTime + /*Real*/(pulse.count) * pulse.period;

  3(1):	pla.coo.con.vol.dynBal.medium.T = pla.coo.con.vol.dynBal.T_start
  5(1):	pla.coo.con.vol.dynBal.medium.Xi[1] = pla.coo.con.vol.dynBal.X_start[1]
  6(2):	pla.coo.con.sta_start.X = pla.coo.con.X_start
  12(1):	pla.hea.eva.vol.dynBal.medium.p = pla.hea.eva.vol.dynBal.p_start
  14(2):	pla.hea.eva.sta_start.X = pla.hea.eva.X_start
  17(1):	vol.dynBal.medium.T = vol.dynBal.T_start
  19(2):	pum.filter.cr = Modelica.Blocks.Continuous.Internal.Filter.base.CriticalDamping(2, pum.filter.normalized)
  20(2):	(pum.filter.r, _, _, _) = Modelica.Blocks.Continuous.Internal.Filter.roots.lowPass(pum.filter.cr, {}, {}, pum.filter.f_cut)
  23(1):	$DER.pum.filter.x[1] = 0.0
  28(1):	$START.pla.coo.port_b1.h_outflow = pla.coo.con.h_outflow_start
  44(1):	pla.sinHea.X_in_internal[2] = pla.sinHea.X[2]
  55(1):	pla.PComHea = pla.QHea_flow / pla.hea.COP
  99(1):	$DER.vol.dynBal.U = vol.dynBal.Hb_flow + heaFlo.Q_flow
  130(1):	pla.hea.port_b2.h_outflow = 1006.0 * pla.hea.eva.vol.dynBal.medium.T_degC * pla.hea.staB2.X[2] + (2501014.5 + 1860.0 * pla.hea.eva.vol.dynBal.medium.T_degC) * pla.hea.eva.vol.dynBal.medium.Xi[1]
  133(1):	pla.QHea_flow = pla.hea.con.outCon.m_flow_pos * pla.hea.con.outCon.dhAct
  135(1):	pla.hea.con.outCon.dhAct = smooth(1, if noEvent(pla.hea.con.outCon.hSet - pla.hea.con.preDro.port_b.h_outflow > pla.hea.con.outCon.deltaH) then pla.hea.con.outCon.hSet - pla.hea.con.preDro.port_b.h_outflow else if noEvent(pla.hea.con.outCon.hSet - pla.hea.con.preDro.port_b.h_outflow < (-pla.hea.con.outCon.deltaH)) then 0.0 else if noEvent(pla.hea.con.outCon.deltaH > 0.0) then 0.25 * (((pla.hea.con.outCon.hSet - pla.hea.con.preDro.port_b.h_outflow) / pla.hea.con.outCon.deltaH) ^ 2.0 - 3.0) * (pla.hea.con.preDro.port_b.h_outflow - pla.hea.con.outCon.hSet) / pla.hea.con.outCon.deltaH * (pla.hea.con.outCon.hSet - pla.hea.con.preDro.port_b.h_outflow) + 0.5 * (pla.hea.con.outCon.hSet - pla.hea.con.preDro.port_b.h_outflow) else 0.5 * (pla.hea.con.outCon.hSet - pla.hea.con.preDro.port_b.h_outflow))
  136(1):	pla.hea.con.outCon.m_flow_non_zero = smooth(1, if noEvent(pum.filter.y - pla.hea.con.outCon.m_flow_small > 0.5 * pla.hea.con.outCon.m_flow_small) then pum.filter.y else if noEvent(pum.filter.y - pla.hea.con.outCon.m_flow_small < (-0.5) * pla.hea.con.outCon.m_flow_small) then pla.hea.con.outCon.m_flow_small else if noEvent(0.5 * pla.hea.con.outCon.m_flow_small > 0.0) then 0.25 * (4.0 * ((pum.filter.y - pla.hea.con.outCon.m_flow_small) / pla.hea.con.outCon.m_flow_small) ^ 2.0 - 3.0) * 2.0 * (pum.filter.y - pla.hea.con.outCon.m_flow_small) / pla.hea.con.outCon.m_flow_small * (pla.hea.con.outCon.m_flow_small - pum.filter.y) + 0.5 * (pum.filter.y + pla.hea.con.outCon.m_flow_small) else 0.5 * (pum.filter.y + pla.hea.con.outCon.m_flow_small))
  175(1):	pla.coo.con.vol.dynBal.U = pla.coo.con.vol.dynBal.m * pla.coo.con.vol.dynBal.medium.u
  177(1):	pla.coo.con.vol.dynBal.mb_flow = 1.184307920059215e-05 * pla.coo.con.vol.dynBal.fluidVolume * $DER.pla.coo.con.vol.dynBal.medium.p
  178(1):	pla.coo.staB1.X[2] = 1.0 - pla.coo.con.vol.dynBal.medium.Xi[1]
  192(1):	pla.coo.sta_a2.T = 273.15 + 0.0002390057361376673 * (if noEvent((-pum.filter.y) > 0.0) then vol.ports[1].h_outflow else pla.port_b.h_outflow)
  193(1):	pla.coo.sta_b1.T = 273.15 + ((if noEvent((-pla.sinCoo.ports[1].m_flow) > 0.0) then pla.sinCoo.ports[1].h_outflow else pla.coo.port_b1.h_outflow) + (-2501014.5) * pla.coo.sta_b1.X[1]) / (1006.0 * (1.0 - pla.coo.sta_b1.X[1]) + 1860.0 * pla.coo.sta_b1.X[1])
  197(1):	pla.coo.sta_a1.X[2] = 1.0 - pla.coo.sta_a1.X[1]
  198(1):	pla.coo.sta_a1.X[1] = if noEvent(pla.coo.m1_flow > 0.0) then pla.sou2.X[1] else 0.01

---- involved set vars (28):
  1:	pum.eff.hydDer[1]:VARIABLE(fixed = false protected = true final = true )  "Coefficients for polynomial of hydraulic efficiency vs. volume flow rate" type: Real [1]
  2:	pum.eff.motDer[1]:VARIABLE(fixed = false protected = true final = true )  "Coefficients for polynomial of motor efficiency vs. volume flow rate" type: Real [1]
  3:	pum.filter.r[2]:VARIABLE(fixed = false protected = true )  type: Real [2]
  4:	pum.filter.r[1]:VARIABLE(fixed = false protected = true )  type: Real [2]
  5:	pum.filter.cr[2]:VARIABLE(fixed = false protected = true )  type: Real [2]
  6:	pum.filter.cr[1]:VARIABLE(fixed = false protected = true )  type: Real [2]
  7:	pla.hea.eva.sta_start.T:VARIABLE(min = 1.0 max = 10000.0 start = 293.15 unit = "K" fixed = false nominal = 300.0 protected = true )  "Temperature of medium" type: Real
  8:	pla.hea.eva.sta_start.X[1]:VARIABLE(min = 0.0 max = 1.0 start = 0.01 unit = "kg/kg" fixed = false nominal = 0.1 protected = true )  "Mass fractions (= (component mass)/total mass  m_i/m)" type: Real [2]
  9:	pla.hea.eva.h_outflow_start:VARIABLE(unit = "J/kg" protected = true )  = 1006.0 * (-273.15 + pla.hea.eva.sta_start.T) * (1.0 - pla.hea.eva.sta_start.X[1]) + (2501014.5 + 1860.0 * (-273.15 + pla.hea.eva.sta_start.T)) * pla.hea.eva.sta_start.X[1]  "Start value for outflowing enthalpy" type: Real
  10:	pla.hea.eva.sta_start.X[2]:VARIABLE(min = 0.0 max = 1.0 start = 0.99 unit = "kg/kg" fixed = false nominal = 0.1 protected = true )  "Mass fractions (= (component mass)/total mass  m_i/m)" type: Real [2]
  11:	pla.coo.con.sta_start.T:VARIABLE(min = 1.0 max = 10000.0 start = 293.15 unit = "K" fixed = false nominal = 300.0 protected = true )  "Temperature of medium" type: Real
  12:	pla.coo.con.sta_start.X[1]:VARIABLE(min = 0.0 max = 1.0 start = 0.01 unit = "kg/kg" fixed = false nominal = 0.1 protected = true )  "Mass fractions (= (component mass)/total mass  m_i/m)" type: Real [2]
  13:	pla.coo.con.h_outflow_start:VARIABLE(unit = "J/kg" protected = true )  = 1006.0 * (-273.15 + pla.coo.con.sta_start.T) * (1.0 - pla.coo.con.sta_start.X[1]) + (2501014.5 + 1860.0 * (-273.15 + pla.coo.con.sta_start.T)) * pla.coo.con.sta_start.X[1]  "Start value for outflowing enthalpy" type: Real
  14:	pla.coo.con.sta_start.X[2]:VARIABLE(min = 0.0 max = 1.0 start = 0.99 unit = "kg/kg" fixed = false nominal = 0.1 protected = true )  "Mass fractions (= (component mass)/total mass  m_i/m)" type: Real [2]
  37:	pum.filter.y:VARIABLE(flow=false min = max(-min(100000.0, 100000.0), max(max(-9.999999999999999e+59, max(-9.999999999999999e+59, max(-100000.0, max(-9.999999999999999e+59, max(-9.999999999999999e+59, max(-100000.0, max(-9.999999999999999e+59, max(-100000.0, max(-9.999999999999999e+59, max(-9.999999999999999e+59, max(-9.999999999999999e+59, max(-9.999999999999999e+59, max(-9.999999999999999e+59, max(-9.999999999999999e+59, max(-9.999999999999999e+59, max(-9.999999999999999e+59, max(-100000.0, max(-100000.0, max(-100000.0, max(-100000.0, max(-100000.0, max(-100000.0, max(-100000.0, max(-100000.0, max(-100000.0, max(-9.999999999999999e+59, -9.999999999999999e+59)))))))))))))))))))))))))), max(max(-100000.0, max(-9.999999999999999e+59, -100000.0)), max(-9.999999999999999e+59, max(-9.999999999999999e+59, -9.999999999999999e+59))))) max = min(-max(-100000.0, -100000.0), min(min(100000.0, min(100000.0, min(100000.0, min(100000.0, min(100000.0, min(100000.0, min(100000.0, min(100000.0, min(100000.0, min(100000.0, min(100000.0, min(100000.0, min(100000.0, min(100000.0, min(100000.0, min(100000.0, min(9.999999999999999e+59, min(9.999999999999999e+59, min(9.999999999999999e+59, min(9.999999999999999e+59, min(9.999999999999999e+59, min(9.999999999999999e+59, min(9.999999999999999e+59, min(100000.0, min(9.999999999999999e+59, min(100000.0, 100000.0)))))))))))))))))))))))))), min(min(100000.0, min(100000.0, 100000.0)), min(100000.0, min(100000.0, 100000.0))))) start = 0.0 unit = "kg/s" nominal = pum.m_flow_nominal protected = true )  "Connector of Real output signal" type: Real
  55:	vol.dynBal.medium.T_degC:VARIABLE(unit = "degC" protected = true )  "Temperature of medium in [degC]" type: Real
  56:	vol.dynBal.medium.T:VARIABLE(min = max(1.0, 1.0) max = min(10000.0, 10000.0) start = 293.15 unit = "K" nominal = 300.0 protected = true )  "Temperature of medium" type: Real
  58:	vol.ports[1].h_outflow:VARIABLE(min = max(max(-10000000000.0, max(-10000000000.0, max(-10000000000.0, -10000000000.0))), max(-10000000000.0, max(-10000000000.0, max(-10000000000.0, max(-10000000000.0, max(-10000000000.0, -10000000000.0)))))) max = min(min(10000000000.0, min(10000000000.0, min(10000000000.0, 10000000000.0))), min(10000000000.0, min(10000000000.0, min(10000000000.0, min(10000000000.0, min(10000000000.0, 10000000000.0)))))) start = pum.h_outflow_start unit = "J/kg" nominal = 1000000.0 )  "Specific thermodynamic enthalpy close to the connection point if m_flow < 0" type: Real [2]
  106:	pla.hea.con.outCon.hSet:VARIABLE(unit = "J/kg" protected = true )  "Set point for enthalpy leaving port_b" type: Real
  109:	pla.hea.con.preDro.dp:VARIABLE(start = 0.0 unit = "Pa" nominal = 30000.0 protected = true )  "Pressure difference between port_a and port_b" type: Real
  110:	pla.hea.con.preDro.port_b.h_outflow:VARIABLE(min = -10000000000.0 max = 10000000000.0 start = 83680.0 unit = "J/kg" nominal = 83680.0 protected = true )  "Specific thermodynamic enthalpy close to the connection point if m_flow < 0" type: Real
  159:	pla.coo.con.vol.dynBal.medium.T_degC:VARIABLE(start = 20.0 unit = "degC" protected = true final = true )  "Temperature of medium in [degC]" type: Real
  163:	pla.coo.con.vol.dynBal.medium.T:VARIABLE(min = max(1.0, 1.0) max = min(10000.0, 10000.0) start = 293.15 unit = "K" nominal = 300.0 protected = true final = true )  "Temperature of medium" type: Real
  166:	pla.coo.con.vol.dynBal.medium.p:DUMMY_STATE(flow=false min = 0.0 start = pla.coo.con.vol.dynBal.p_start unit = "Pa" nominal = 100000.0 stateSelect=StateSelect.never protected = true )  "Absolute pressure of medium" type: Real
  185:	pla.coo.sta_a1.X[1]:VARIABLE(min = 0.0 max = 1.0 start = 0.01 unit = "kg/kg" nominal = 0.1 protected = true )  "Mass fractions (= (component mass)/total mass  m_i/m)" type: Real [2]
  189:	$START.pla.coo.port_b1.h_outflow:VARIABLE(min = max(-10000000000.0, max(-10000000000.0, max(-10000000000.0, max(-10000000000.0, max(-10000000000.0, -10000000000.0))))) max = min(10000000000.0, min(10000000000.0, min(10000000000.0, min(10000000000.0, min(10000000000.0, 10000000000.0))))) unit = "J/kg" fixed = false nominal = 45300.945 protected = true )  "Specific thermodynamic enthalpy close to the connection point if m_flow < 0" type: Real
  190:	pla.coo.port_b1.h_outflow:VARIABLE(min = max(-10000000000.0, max(-10000000000.0, max(-10000000000.0, max(-10000000000.0, max(-10000000000.0, -10000000000.0))))) max = min(10000000000.0, min(10000000000.0, min(10000000000.0, min(10000000000.0, min(10000000000.0, 10000000000.0))))) start = pla.coo.con.h_outflow_start unit = "J/kg" nominal = 45300.945 protected = true )  "Specific thermodynamic enthalpy close to the connection point if m_flow < 0" type: Real
  198:	pla.sta_b.p:VARIABLE(min = max(max(0.0, 0.0), max(max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, 0.0))))))))))), max(max(0.0, 0.0), max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, max(0.0, 0.0))))))))))))))))))) max = min(min(100000000.0, 100000000.0), min(min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, 100000000.0))))))))))), min(min(100000000.0, 100000000.0), min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, min(100000000.0, 100000000.0))))))))))))))))))) start = 300000.0 unit = "Pa" nominal = 100000.0 )  "Absolute pressure of medium" type: Real
--------------------------------------------------

[Karim Adbdelhak]

To understand this output i should add that there are 26 scalarized equations counted. The number of shown equations (compact as array equations) is 25. For each individual equation we see the equation index and in brackets the size of the equation. As we can see there are multiple two dimensional equations but only one extra scalarized equation, this means that some of these array equations are only half part of this set and the other half seemed to be perfectly matched in the system (i don't know how to interpret this, but that is how to read this).

[Francesco Casella]

@Karim, I had a look at this output, but I am a bit confused about what it means.
The following variables do not show up in the equations:

pum.eff.hydDer[1]
pum.eff.motDer[1]

the following ones are used in binding equations but seem not to have any equation to define them

pla.hea.eva.sta_start.T
pla.hea.eva.sta_start.X[1]
pla.hea.eva.sta_start.X[2]
pla.coo.con.sta_start.T
pla.coo.con.sta_start.X[1]
pla.coo.con.sta_start.X[2]

Last, but not least, equation

pla.PComHea = pla.QHea_flow / pla.hea.COP

contains three variables, none of which show up in the list of variables.

I guess something went badly wrong here. Do I miss something?

[Francesco Casella]

Ticket retargeted after milestone closed

[Francesco Casella]

Smaller example giving the same error: ​Buildings.Fluid.HeatExchangers.DXCoils.AirCooled.Examples.SingleSpeed

Error: The initialization problem of given system is mixed-determined.
It is under- as well as overdetermined and the mixed-determination-index is too high. [index > 10]

[Francesco Casella]

The two test cases reported here work after merging ​PR #7903:

• ​Buildings.Fluid.HeatExchangers.DXCoils.AirCooled.Examples.SingleSpeed
• ​Buildings.Experimental.DistrictHeatingCooling.Plants.Validation.Plant_Carnot_T_ClosedLoop