Check models that fail as actuator model has been updated

[mwetter]

The new actuator model (#4000) causes some models to fail to simulate, see conf.yml in
https://github.com/lbl-srg/modelica-buildings/pull/4000/files#diff-27f3a54beb1453b97cf6c7dfd9fa86b042c9c7e332cd38ca55ea859c936120a8

This issue is to see if there is something that need to be improved in these models.

[mwetter]

@AntoineGautier : Could you please have a look into why these models fail with the new actuator model. These are from the HP template and the hydronic configurations. I believe the hydronic configuration have a pump that works against a fully closed valve.
I am not sure why Buildings.Templates.Plants.HeatPumps.Components.Validation.ValvesIsolation fails with OCT, but at the time of when OCT stalls, there is a discontinuity in slope:

[AntoineGautier]

Models that fail to simulate with the change

- model_name: Buildings.DHC.Plants.Combined.Validation.AllElectricCWStorage
  optimica:
    comment:  "Simulation times out"
    simulate: false
- model_name: Buildings.Fluid.HydronicConfigurations.ActiveNetworks.Examples.InjectionTwoWayCheckValve
  openmodelica:
    comment: 'Singular system of equations when both load valves are fully closed but pump is running'
    simulate: false
- model_name: Buildings.Templates.Plants.HeatPumps.Components.Validation.ValvesIsolation
  optimica:
    comment: Possible chattering detected at t = 3.120112e+03 in state event
    simulate: false

Simulating these models with Dymola also yields NL solver failures, but the solver can recover from them and the simulation ultimately succeeds.

We can already notice that all these models include one or more instances of Buildings.Fluid.FixedResistances.CheckValve and that the NL solver issues disappear when these instances are removed, or even replaced by a fixed resistance.

DHC.Plants.Combined.Validation.AllElectricCWStorage

Locally, the simulation succeeds with OCT v1.51.6 (same as in .travis.yml) on Ubuntu 22.04 with the attached Python dependencies requirements_1.51.6.txt that match the requirements from "OPTIMICA Compiler Toolkit Version 1.51.6" Appendix C.

Tolerances (relative)    : 1e-06
Simulation interval    : 0.0 - 345600.0 seconds.
Elapsed simulation time: 70.70019977999982 seconds.

@mwetter Can you check on your system if the simulation succeeds with OCT v1.51.6?

Fluid.HydronicConfigurations.ActiveNetworks.Examples.InjectionTwoWayCheckValve

Locally, the simulation stalls with OMC (v1.24.0~dev.beta.0-5-g904c4c7 VS v1.23.0-1 in .travis.yml) at around 6E4 s when loa1.yVal_actual ~ 0.63 and con.yVal_actual ~ 0.59 — so not in a situation where the pump work against closed valves. Anyway both primary and secondary loops have an end-of-line bypass that prevents the pumps from being dead-headed.

The time of failure matches a local minimum of the check valve ∆p:

With the following change in Buildings.Fluid.FixedResistances.CheckValve, the simulation succeeds with OMC.
This is hard to understand though, as a^3*(10+a*(-15+6*a)) is a monotonous C∞ function that equals 0 at a=0 and 1 at a=1 and f"(0) = f"(1) = 0 so there is continuity of the second derivative with the min and max bounding functions.

  // cv = smooth(2, max(0, min(1, a^3*(10+a*(-15+6*a)))));
  cv = Buildings.Utilities.Math.Functions.smoothLimit(a^3*(10+a*(-15+6*a)), 0, 1, 1E-4);

### STATISTICS ###
timer
    0.112537s          reading init.xml
   0.0179286s          reading info.xml
  0.00424254s [  1.0%] pre-initialization
    0.093949s [ 21.3%] initialization
 0.000161417s [  0.0%] steps
   0.0205117s [  4.7%] solver (excl. callbacks)
  0.00305929s [  0.7%] creating output-file
   0.0025535s [  0.6%] event-handling
  0.00428121s [  1.0%] overhead
     0.31148s [ 70.8%] simulation
    0.440238s [100.0%] total

Templates.Plants.HeatPumps.Components.Validation.ValvesIsolation

Once again, the check valve component appears to cause the NL solver convergence failure with OCT.
The NL equation block 24 is reported in the log, which contains:

cheValChiWat[2].m_flow	             
cheValChiWat[1].m_flow	             
cheValChiWat[1].cv	                 
cheValHeaWat[2].dp	                 
valIsoSep.valChiWatHpInlIso[2].m_flow
cheValChiWat[1].port_a.p	         
cheValHeaWat[1].dp

This calls for improving the check valve model. We could:

1. Add an option linearized to allow using a linear relationship between dp and m_flow. This is a generic solution that makes the simulation succeeds with OCT for this particular model, and that reduces the simulation time by 30 % for DHC.Plants.Combined.Validation.AllElectricCWStorage with OCT. The error on the pump power due to approximating the check valve pressure drop is obviously marginal on large fluid systems.

2. Add an option from_dp to allow solving for dp instead of m_flow (using basicFlowFunction_m_flow instead of basicFlowFunction_dp). This makes the simulation succeeds with OCT for this particular model.
   However, this is not a generic solution as, with this change, other models fail to simulate with Dymola, such as DHC.Plants.Combined.Validation.AllElectricCWStorage.

3. Use an "ideal diode" approach such as in Modelica.Fluid.Machines.BaseClasses.PartialPump — although I don't know the robustness of this model when used in complex fluid systems.
   Note that this model also scales the flow characteristic to represent multiple pumps in parallel, which seems a good alternative to connecting multiple instances of Buildings.Fluid.Movers.SpeedControlled_y and Buildings.Fluid.FixedResistances.CheckValve as we currently do.

@mwetter Let me know what you think of these options.

[mwetter]

@AntoineGautier The model "DHC.Plants.Combined.Validation.AllElectricCWStorage" works on my installation.

I would think 1 is the most promising option, followed by 2. I don't see why 3 should be favorable to what we have. Another option would be to make the regularization near 0 less abrupt.

@justnielsen: Are you using the check valve of Buildings in conjunction with the actuator that has linear travel time, or do you have an alternate implementation of the check valve with which you have good numerical performance?

[justnielsen]

@mwetter we have just recently updated our valve and pump models to reflect the actuator changes in Buildings Library. For the majority of the many models we tested, the conversion went without problems. But a few of our very large models (100–200,000 equations) needed some tweaking. During this we did not consider looking at the check valves but generally solved the problems (failing to converge/initialize + slow simulation) by adding/removing dynamics in T-junctions (Buildings.Fluid.FixedResistances.Junction).

Regarding our check valve model it is an unmodified extension of (Buildings.Fluid.FixedResistances.CheckValve). We have used it with great success in many models — typically downstream of pumps — and have (AFAIK) not had any problems running pumps against closed valves downstream of the check valves. This would be the case during start-up of a pump/valve controlled system where the pump must run at minimum speed before the valve can open.

However, we often experience nonlinear system convergence errors in Dymola but unfortunately, our analyses of the numerical behaviour of these models have never been near as qualified and rigorous as the one you present in this thread.

[AntoineGautier]

@mwetter I realized that even after substituting the call to basicFlowFunction_dp with a linear relationship between m_flow and dp, we still have a 5-degree polynomial relationship between the flow coefficient and dp that is chosen to ensure a continuous second derivative with the bounding flow coefficients.
Overall, I would favor using a quintic Hermite spline interpolation (C2) between the closed valve flow function and the fully open valve flow function — similar approach as in Fluid.Actuators.Dampers.PressureIndependent. I quickly implemented this approach and it eliminates all convergence issues on all the models from this thread (with minor but favorable impact on simulation time). Even recoverable NL solver failures with Dymola are eliminated.

If you agree, I can submit a PR to incorporate this reformulation into the IBPSA library.

[mwetter]

@AntoineGautier : Using a spline is probably better here. Perhaps the solver extrapolated on the polynomial while trying to find a solution, and polynomials can rapidly have huge slopes when used outside their intended range.

Please do the PR.