# 3. Simulation Performance¶

This section provides tips for how to implement numerically efficient models and for how to debug models.

All Modelica simulation environments have debugging and profiling capabilities that give insight into where computing time is spent and how to speed up simulations. The use of these features varies among the different Modelica tools. We recommend users read the documentation of their respective simulation environment, take opportunity of training that is offered on request by many tool providers, and contact their support for tool related questions.

## 3.1. Unstable control loops¶

Most solvers in Modelica tools adapt the time step during the simulation in order to accurately resolve the dynamics of the model. If feedback control loops are unstable, then trajectories can have high frequency oscillations, just as in a real HVAC system. The solver is then required to simulate with small time steps, which increases the simulation time.

In large models, to detect which control loop is unstable, it can be helpful to log which state variables dominate the integration error or the integrator time step control. This usually points to the control loop that is unstable. There are various methods for tuning a controller, and the documentation of Buildings.Controls.OBC.CDL.Reals.PIDWithReset outlines one approach for tuning the gains of a PI-controller.

## 3.2. Control input signal of equipment¶

Most equipment models that take a real-valued control signals as an input, such as
flow machines (fans and pumps), have a boolean parameter
`filteredSpeed`

, and all actuators have a boolean parameter
`filteredOpening`

.
If set to `true`

, which is the default setting, then the control input signal is sent to
a 2nd order low pass filter that changes a step signal to a smooth signal.
This typically improves the robustness of the simulation.

To see the effect of the filter, consider the model below
in which `fanS`

is configured with
`filteredSpeed=false`

, and `fanC`

is configured with
`filteredSpeed=true`

.
Both fans are connected to a step input signal.
The configuration of `fanS`

causes the fan speed to instantly change from 0 to 1. In large system models, this can lead to high computing time or to convergence problems. The `fanC`

avoids this problem because the speed of the fan varies continuously, thereby making it easier for the solver to compute a solution. In this model, we set the parameter
`raiseTime=30`

seconds.

For fans and pumps, the dynamics introduced by the filter can be thought of as approximating
the rotational inertia of the fan rotor and the inertia of the fluid in the duct or piping network.
The default value is `raiseTime=30`

seconds.

For actuators, the raise time approximates the travel time of the valve lift.
The default value is `raiseTime=120`

seconds.

Note

When changing `filteredSpeed`

(or `filteredOpening`

),
or when changing the value of `raiseTime`

, the dynamic
response of the closed loop control changes. Therefore,
control gains may need to be retuned to ensure satisfactory
closed loop control performance.

For further information, see the User’s Guide of the flow machine package, and the User’s Guide of the actuator package.

## 3.3. Fluid flow systems¶

### 3.3.1. Breaking algebraic loops¶

In fluid flow systems, flow junctions where mass flow rates separate and mix can couple non-linear systems of equations.
This leads to larger systems of coupled equations that need to be solved,
which often causes larger computing time and can sometimes cause convergence problems.
To decouple these systems of equations, in the model of a flow junction
(Buildings.Fluid.FixedResistances.Junction),
or in models for fans or pumps (such as the model
Buildings.Fluid.Movers.SpeedControlled_y),
the parameter `dynamicBalance`

can be set to `true`

.
This adds a control volume at the fluid junction that can decouple the system of equations.

### 3.3.2. Reducing nonlinear equations of serially connected flow resistances¶

In fluid flow systems, if multiple components are connected in series, then computing the pressure drop due to flow friction in the individual components can lead to coupled nonlinear systems of equations. While this is no problem for small models, the iterative solution can lead to higher computing time, particularly in large models where other equations may be part of the residual function.

For illustration, consider the simple system shown below in which the flow resistances `res1`

and `res2`

compute the mass flow rate as
\(\dot m = k \, \sqrt{\Delta p}\) if the parameter `from_dp`

is set to `true`

, or otherwise compute the pressure drop between their inlet and outlet as \(\Delta p = (\dot m / k)^2\). (Both formulations are implemented using regularization near zero.)

Depending on the configuration of the individual component models, simulating this system model may require the iterative solution of a nonlinear equation to compute the mass flow rate or the pressure drop. To avoid a nonlinear equation, use any of the measures below.

Set the parameter

`res2(dp_nominal=0)`

, and add the pressure drop to the parameter`dp_nominal`

of the model`res1`

. This will eliminate the equation that computes the flow friction in`res2`

, thereby avoiding a nonlinear equation. The same applies if there are multiple components in series, such as a pre-heat coil, a heating coil and a cooling coil.Set

`from_dp=false`

in all components, which is the default setting. This will cause Modelica to use a function that computes the pressure drop as a function of the mass flow rate. Therefore, a code translator is likely to generate an equation that solves for the mass flow rate, and it then uses the mass flow rate to compute the pressure drop of the components that are connected in series.

Control valves also allow lumping the pressure drop into the model of the valve. Consider the situation where a fixed flow resistance is in series with a control valve as shown below.

Suppose the parameters are

```
Buildings.Fluid.FixedResistances.PressureDrop res(
redeclare package Medium = Medium,
m_flow_nominal=0.2,
dp_nominal=10000);
Buildings.Fluid.Actuators.Valves.TwoWayLinear val(
redeclare package Medium = Medium,
m_flow_nominal=0.2,
dpValve_nominal=5000);
```

To avoid a nonlinear equation, the flow resistance could be deleted as shown below.

If the valve is configured as

```
Buildings.Fluid.Actuators.Valves.TwoWayLinear val(
redeclare package Medium = Medium,
m_flow_nominal=0.2,
dpValve_nominal=5000,
dpFixed_nominal=10000);
```

then the valve will compute the composite flow coefficient \(\bar k\) as

where \(k_v(y) = \dot m(y)/\sqrt{\Delta p}\) is the flow coefficient of the valve at the lift \(y\), and
\(k_f\) is equal to the ratio `m_flow_nominal/sqrt(dpFixed_nominal)`

.
The valve model then computes the pressure drop using \(\bar k\) and the same equations as described above for the fixed resistances.
Thus, the composite model has the same valve authority and mass flow rate, but a nonlinear equation can be avoided.

For more details, see the User’s Guide of the actuator package.

### 3.3.3. Prescribed mass flow rate¶

For some system models, the mass flow rate can be prescribed by using an idealized pump or fan (model Buildings.Fluid.Movers.FlowControlled_m_flow) or a source element that outputs the required mass flow rate (such as the model Buildings.Fluid.Sources.MassFlowSource_T). Using these models avoids having to compute the intersection of the fan curve and the flow resistance. In some situations, this can lead to faster and more robust simulation.

Warning

If you prescribe the mass flow rate, make sure the pump (or fan) does not work agains a closed valve (or damper). Otherwise, it will force the flow through the component, which leads to very large pressure drops.

### 3.3.4. Avoiding events¶

In Modelica, the time integration is halted whenever a Real elementary
operation such as \(x>y\), where \(x\) and \(y\) are variables of type `Real`

,
changes its value. In this situation,
an event occurs and the solver determines a small interval in time in which
the relation changes its value. This can increase computing time.
An example where such an event occurs is the following relation
that computes the enthalpy of the medium that streams through `port_a`

as

```
if port_a.m_flow > 0 then
h_a = inStream(port_a.h_outflow);
else
h_a = port_a.h_outflow;
end if;
```

or, equivalently,

```
h_a = if port_a.m_flow > 0 then inStream(port_a.h_outflow) else port_a.h_outflow;
```

When simulating a model that contains such code, a time integrator
will iterate to find the time instant where `port_a.m_flow`

crosses zero.
If the modeling assumptions allow approximating this equation in
a neighborhood around `port_a.m_flow=0`

, then replacing this equation
with an approximation that does not require an event iteration can
reduce computing time. For example, the above equation could be
approximated as

```
T_a = Modelica.Fluid.Utilities.regStep(
port_a.m_flow, inStream(port_a.h_outflow), port_a.h_outflow,
m_flow_nominal*1E-4);
```

where `m_flow_nominal`

is a parameter that is set to a value that
is close to the mass flow rate that the model has at full load.
If the magnitude of the flow rate is larger than 1E-4 times the
typical flow rate, the approximate equation is the same as the exact equation,
and below that value, an approximation is used. However, for such small
flow rates, not much energy is transported and hence the error introduced
by the approximation is generally negligible.

In some cases, adding dynamics to the model can further improve
the computing time, because the return value of the function
Modelica.Fluid.Utilities.regStep()
above can change abruptly if its argument `port_a.m_flow`

oscillates in the range of
`+/- 1E-4*m_flow_nominal`

,
for example due to numerical noise.
Adding dynamics may be achieved using a formulation such as

```
hMed = Modelica.Fluid.Utilities.regStep(
port_a.m_flow, inStream(port_a.h_outflow), port_a.h_outflow,
m_flow_nominal*1E-4);
der(h)=(hMed-h)/tau;
```

where `tau`

>0 is a time constant. See, for example,
Buildings.Fluid.Sensors.SpecificEnthalpyTwoPort
for a robust implementation.

Note

In the package
Buildings.Utilities.Math
the functions and blocks whose names start with `smooth`

can be used to avoid events.

## 3.4. Examples for how to debug and correct slow simulations¶

### 3.4.1. State events¶

This section shows how a simulation that stalls due to events can be debugged
to find the root cause, and then corrected.
While the details may differ from one tool to another, the principle is the same.
In our situation, we attempted to simulate `Buildings.Examples.DualFanDualDuct`

for one year in Dymola 2016 FD01 using the model from Buildings version 3.0.0.
We run

```
simulateModel("Buildings.Examples.DualFanDualDuct.ClosedLoop",
stopTime=31536000, method="radau",
tolerance=1e-06, resultFile="DualFanDualDuctClosedLoop");
```

and plotted the computing time and the number of events. Around \(t=0.95e7\) seconds, there was a spike as shown in the figure below.

As the number of events increased drastically, we enabled in Dymola in
Simulation -> Setup, under the tab Debug the entry Events during simulation
and simulated the model from
\(t=0.9e7\) to \(t=1.0e7\) seconds. It turned out that setting the start time
to \(t=0.9e7\) seconds was sufficient to reproduce the behavior;
otherwise we would
have had to set it to an earlier time.
Inspecting Dymola’s log file `dslog.txt`

when the simulation stalls shows that its last entries
are

```
Expression TRet.T > amb.x_pTphi.T became true ( (TRet.T)-(amb.x_pTphi.T) = 2.9441e-08 )
Iterating to find consistent restart conditions.
during event at Time : 9267949.854873843
Expression TRet.T > amb.x_pTphi.T became false ( (TRet.T)-(amb.x_pTphi.T) = -2.94411e-08 )
Iterating to find consistent restart conditions.
during event at Time : 9267949.855016639
Expression TRet.T > amb.x_pTphi.T became true ( (TRet.T)-(amb.x_pTphi.T) = 2.94407e-08 )
Iterating to find consistent restart conditions.
during event at Time : 9267949.855208419
Expression TRet.T > amb.x_pTphi.T became false ( (TRet.T)-(amb.x_pTphi.T) = -2.94406e-08 )
Iterating to find consistent restart conditions.
during event at Time : 9267949.855351238
```

Hence, there is an event every few milliseconds, which explains
why the simulation does not appear to be progessing.
The solver does the right thing, it stops
the integration, handles the event, and restarts the integration, just to encounter
another event a few milliseconds later.
Hence, we go back to our system model and
follow the output signal of `TRet.T`

of the
return air temperature sensor,
which shows that it is used in the economizer control
to switch the sign of the control gain because the economizer can provide heating or cooling,
depending on the ambient and return air temperature. The problematic model is
shown in the figure below.

The events are triggered by the inequality block which changes the control, which then in turn seems to cause a slight change in the return air temperature, possibly due to numerical noise or maybe because the return fan may change its operating point as the dampers are adjusted, and hence change the heat added to the medium. Regardless, this is a bad implementation that also would cause oscillatory behavior in a real system if the sensor signal had measurement noise. Therefore, this equality comparison must be replaced by a block with hysteresis, which we did as shown in the figure below. We selected a hysteresis of \(0.2\) Kelvin, and now the model runs fine for the whole year.

### 3.4.2. State variables that dominate the error control¶

In a development version of the model
`Buildings.Examples.DualFanDualDuct.ClosedLoop`

(commit ef410ee),
the simulation time was very slow during part of the
simulation, as shown in Fig. 3.9.

The number of state events did not increase in that time interval. To isolate the problem, we enabled in Dymola under Simulation -> Setup the option to log which states dominate the error (see Debug tab).

Running the simulation again gave the following output:

```
Integration terminated successfully at T = 1.66e+07
Limit stepsize, Dominate error, Exceeds 10% of error Component (#number)
0 1 6 cooCoi.temSen_1.T (# 1)
36 0 140 cooCoi.temSen_2.T (# 2)
37 0 0 cooCoi.ele[1].mas.T (# 3)
45 0 0 cooCoi.ele[2].mas.T (# 4)
51 0 0 cooCoi.ele[3].mas.T (# 5)
53 0 0 cooCoi.ele[4].mas.T (# 6)
13555 13201 19064 fanSupHot.filter.x[1] (# 7)
11905 2170 12394 fanSupHot.filter.x[2] (# 8)
400 47 419 fanSupCol.filter.x[1] (# 9)
420 71 521 fanSupCol.filter.x[2] (# 10)
5082 2736 6732 fanRet.filter.x[1] (# 11)
1979 25 4974 fanRet.filter.x[2] (# 12)
38 0 3 TPreHeaCoi.T (# 13)
30 0 1 TRet.T (# 14)
38 0 3 TMix.T (# 15)
80 0 0 TCoiCoo.T (# 16)
305 22 275 cor.vavHot.filter.x[1] (# 18)
```

Hence, the state variables in the highlighted lines
limit the step size significantly more often than other variables.
Therefore, we removed these state variables
by setting in the fan models the parameter `filteredSpeed=false`

.
After this change, the model simulates without problems.

## 3.5. Numerical solvers¶

Dymola 2021 is configured to use dassl as a default solver with a tolerance of 1E-4. We recommend to change this setting to radau with a tolerance of around 1E-6, as this generally leads to faster and more robust simulation for thermo-fluid flow systems.

Note that this is the error tolerance of the local integration time step. Most ordinary differential equation solvers only control the local integration error and not the global integration error. As a rule of thumb, the global integration error is one order of magnitude larger than the local integration error. However, the actual magnitude of the global integration error depends on the stability of the differential equation. As an extreme case, if a system is chaotic and uncontrolled, then the global integration error will grow rapidly.