Implement linearization algorithm for damper and valve models

[AntoineGautier]

Difficulties are encountered by users when it comes to properly sizing dampers and valves. This results in poor authorities, unstable control loops and eventually in various strategies to deal with the non linearity of the controlled system e.g. in this post. Another typical case might be the control of the secondary fluid outlet temperature of an heat exchanger by actuating the valve position on the primary side.

The development aims at providing a built-in option in the actuator models to linearize the controlled system in its whole. It only addresses static non linearities (deviation of the whole system characteristic from linear trend) as opposed to dynamic non linearities (related to actuators' dynamics).

[AntoineGautier]

@Mathadon @mwetter Hi Filip Related to the current issue we ran some tests with Buildings.Fluid.Actuators.Dampers.PressureIndependent. See in the forked branch here: issue1298_tests The test consists in modeling a pressure independent damper (dp_nominal = 10 Pa) with an upstream fixed flow resistance (dp_nominal = 90 Pa) and a fan operating at constant speed (dp_nominal = 100 Pa). At low flow rate it shows a significant discrepancy between setpoint and actual flow rate. This happens because the pressure drop is driven by the fan and largely exceeds the pressure drop of the fully open damper (dp_min in the code) thus leading to m_flow_set + (dp - dp_min) * l2 * m_flow_nominal / dp_nominal significantly higher than m_flow_set. However it seems that such a case would be typical of VAV systems with pressure independent terminal box (system pressure drop much higher than terminal branch pressure drop). If we set l2 to a very low value it lessens the discrepancy but then the nominal operating point is impacted with option from_dp = false. A simple model has been tested in comparison: Buildings.Fluid.Actuators.Dampers.VAVBoxExpLinTest. It presents a good behavior but to my opinion there is no way to implement a from_dp = false option for this algorithm. Indeed, whatever the pressure drop, the damper will adjust and generate the flow resistance needed to reach the setpoint flow rate (within the limit of its fully open characteristic). Could we get in touch and discuss further this test case?

[Mathadon]

@AntoineGautier thanks for reporting this.

From the model equations:

l2=0.01
coeff1 = l2/dp_nominal*m_flow_nominal
m_flow_y2 = m_flow_set + coeff1*max(dp_x,dp_x2);

it would seem that your flow rate will be m_flow_set + 0.01*100/10*m_flow_nominal=1.1*m_flow_nominal for a pressure difference that is 10 times the nominal pressured rop. I don't consider this to be unrealistic? Or are you seeing different results? As you explain, this difference can be reduced by using the parameter l2. If that does not work for from_dp is false, can you include a figure that illustrates more precisely the case where the results are unexpected? Then I'm sure we're talking about the same thing.

Thanks!

[AntoineGautier]

@Mathadon I agree with the order of magnitude. To me the default settings should indeed handle pressure drops much higher than dp_min (pressure drop of fully open damper + fixed resistance) with the damper model still being able to provide the setpoint flow rate. If we consider a VAV system with constant static pressure control by the supply fan and pressure drop nominal values of: 800 Pa for the main branch and 150 Pa for the VAV box. At 50% nominal flow rate on the main branch and 30% nominal flowrate on the VAV box we would have dp ~ 50 * dp_min. Here under is the simulation plot for (see Test):

• Damper control signal ramp from 0 to 1 in 1 s.
• Nominal: system pressure drop 100 Pa, damper pressure drop 10 Pa, flow rate 1 kg/s.
• Case damPreIndDpFan: default values for advanced parameters (l2 = 0.01 delta_x = 0.02 from_dp = true), system pressure drop provided by a fan at constant speed (so following fan flow characteristic i.e. decreasing with damper opening).
• Case damPreIndDpCst: same as before except for the pressure drop being constant (100 Pa).
• Case damPreIndDpCstl2: same as before except parameter l2 = 1e-6: provides the expected results.
• Case damPreIndDpCstl2FromDp: same as before except parameter from_dp = false: flow rate bias at nominal conditions.
• The bias at nominal conditions from the last case can eventually be lessened with delta_x = 1e-3. The results are then similar to damPreIndDpCstl2.

So what do you think of:

1. Updating the models by modifying default values of l2 and delta_x. Then we should probably perform some additional tests to cover a larger range of potential operating points.
2. Switching to a simpler model like this one: VAVBoxExpLin. Which has the advantage of providing an "exact" solution but lacks an implementation of from_dp = false option.

[mwetter]

I am concerned about setting l2 = 1e-6 as this leakage flow is in the same, or even smaller, order as the typically selected solver tolerance of 1E-5 to 1E-6.

[Mathadon]

I'm not sure whether using 1e-6 would be a problem to be honest. Its relevance is mostly with respect to solving algebraic loops? Assuming from_dp=True and the case where the equations ends up in an algebraic loop in this regularization regime we have

m_flow = m_flow_set + coeff1*(dp-dp_min);

according to wikipedia the condition number of a single equation f(x) is x f'(x)/f(x) which is for dp=1000 Pa, dp_min=500Pa and m_flow_set=0.1kg/s : 1000*coeff1/f(x) with f(x)~0.1 since we chose l2 sufficiently small. Assuming coeff1~1e-6*0.1/500 we get a condition number of 1e-6. This causes us to loose about 6 digits of precision in the Newton solver (I think?), which leaves us 10 digits for reaching an accuracy of about 8 digits for the typical Newton solver? Sounds good to me!

If you agree with this analysis I propose to run unit tests with l2=1e-5 or 1e-6 and see what happens. Note that the valve leakage factor of a regular valve is already much smaller than 1e-2.

[AntoineGautier]

We are not sure how the condition number considerations hold in case of nested iteration schemes. On top of unit tests it would be great if you had the opportunity to run tests on large whole system models where both accuracy resulting from the new parameters values and convergence issues will be at stake.

[AntoineGautier]

@Mathadon We came up with an updated version of the PresssureIndependent model using quinticHermite function for smoothing the transitions between the different flow regimes, see: https://github.com/AntoineGautier/modelica-buildings/blob/issue1298_tests/Buildings/Fluid/Actuators/Dampers/VAVBoxPreIndTest.mo. The model is further described and tested here: https://antoinegautier.github.io/Notes/linear-actuators.html#pressure-independent. On this page you might also have a look at our proposal regarding the updates to the Buildings library. The new version offers the advantage of not relying on computational parameters (l2 and delta_x) any more and improving the leakage modeling. What is your opinion on this?

[Mathadon]

@AntoineGautier I had a quick look at it, looks clean!

It looks like you simply removed l2 and delta_x? There thus exists a region where m_flow=y_actual * m_flow_nominal. As such we have d m_flow / d dp = 0. If this equation ends up in an algebraic loop where m_flow is prescribed (e.g. to m_flow_prescribed using IDEAS.Fluid.Sources.MassFlowSource_T) and dp is an iteration variable, the Newton solver will not be able to find a descent direction for dp that reduces the residual equations m_flow - m_flow_prescribed. This was the initial rationale for introducing l2 and deltax. I'm not sure whether this would often cause problems in practice though.. It's perhaps worth making the example that I describe and to see what happens.

[mwetter]

@Mathadon : That is a good catch! Yes, we should have a dependency so that d m_flow / d dp has a non-zero derivative over the whole operating region.

[AntoineGautier]

@mwetter I have just pushed the branch issue1298_linearActuators with the updates on the dampers models discussed in Aachen:

• VAVBoxExponential is retired. The option for the fixed resistance has been integrated into Exponential. Former instances of Exponential must be updated with Exponential(dp_nominal=0, dp_nominalIncludesDamper=false, ...). There is a Python script which automates the required conversions on the whole library: should it be included into the library as well?
• The option to linearize the flow characteristics of the damper plus fixed resistance has been added.
• New limiting values for parameters k0 and k1 are used in PartialDamperExponential.
• PressureIndependent: the fractional opening is now connected to the y_actual output for consistency with other dampers models. However I could not find a proper way to avoid the steep transition of the opening from 1 to 0 when y=0 and dp transitions from negative to positive values. Indeed the actuator's dynamics is represented in a simplified way using the filtered input control value to compute the required mass flow rate y_filtered * m_flow_nominal. There is one main drawback with this approach: at fixed control input values and varying pressure drop the damper reacts instantaneously (the mass flow rate remains constant meaning the opening speed is not limited). So with the current implementation we cannot rely on the opening dynamics to handle the y=0 issue. I have run some preliminary tests with an alternative implementation where the required opening is first computed (with basicFlowFunction_inv(m_flow=y * m_flow_nominal...) and exponentialDamper_inv) then filtered with ActuatorSignal and eventually the flow calculation is performed with Exponential model. However the computational performance is significantly degraded. I am also reluctant to add an additional condition to force a closed position if the control signal is below a threshold value because it just seems like shifting the issue. I would need to meet with you and discuss the best option.