Weather data and solar radiation¶
Thermal circuit and state-space representation for a thermal circuit with capacities in every node: simple wall¶
Thermal circuit and state-space representation for a thermal circuit with capacities in some nodes: cubic building¶
Thermal circuits assembling¶
Switch between models: heating & cooling and free-running¶
Objectives:
- Implement controller with deadband.
- Switch between modeles during the numerical integration.
Figure 1. Thermal circuit for the cubic building.
In this example, we will consider that the controller, modeled by the conductance $G_{11}$ and the temperature source $T_{i,sp}$, switches the value of conductance $G_{11}$ between $G_{11} = 0$, i.e. free-floating, and $G_{11} \rightarrow \infty$, i.e. perfect controller (Figure 1). Therefore, two models will be constructed and the algorithm will switch between these two models during the numerical integration:
- if $\theta_i < T_{i,sp}$ or $\theta_i > T_{i,sp} + \Delta \theta$, then the controller is in action, $G_{11} \equiv K_p \rightarrow \infty$;
- else, the controller is shut off, $G_{11} \equiv K_p = 0$,
where:
- $\theta_i$ is the indoor temperature, °C;
- $T_{i,sp}$ - indoor temperature setpoint, °C;
- $\Delta \theta$ - indoor temperature deadband, °C;
- $K_p$ - static gain of the proportional controller, W/K.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import dm4bem
import tuto
Note: The module tuto.py contains the implementation of the model as presented in this notebook. If the model in this notebook is changed, the same changes need to be done in tuto.py.
filename = './weather_data/FRA_Lyon.074810_IWEC.epw'
start_date = '2000-05-01 00:00:00'
end_date = '2000-05-22 00:00:00'
We chose the [setpoint](https://en.m.wikipedia.org/wiki/Setpoint_(control_system) and the dead-band for the indoor temperature.
Tisp = 20 # indoor setpoint temperature, °C
Δθ = 5 # temperature deadband, °C
Model for building in free-running¶
The building is in free running (the indoor temperature is free-floating) if the P-controller is not effective, i.e. $K_p \rightarrow 0$.
Kpf = 1e-3 # no controller Kp -> 0
TCa = tuto.thermal_circuit(Kpf)
[Af, Bf, Cf, Df] = dm4bem.tc2ss(
TCa['A'], TCa['G'], TCa['C'], TCa['b'], TCa['f'], TCa['y'])
dt_max = min(-2. / np.linalg.eig(Af)[0])
print(f'Maximum time step in free-floating is {dt_max:.2f} s')
Maximum time step in free-floating is 479.08 s
def P_control(filename, start_date, end_date, dt,
As, Bs, Cs, Ds, Kp, Tisp):
t, u, data = tuto.inputs(filename, start_date, end_date, dt, Tisp)
# Initialize temperature vector
θ = np.full([As.shape[0], u.shape[0]], np.inf)
θ[:, 0] = Tisp * np.ones(As.shape[0])
I = np.eye(As.shape[0])
# Simulation in time
for k in range(u.shape[0] - 1):
θ[:, k + 1] = (I + dt * As) @ θ[:, k]\
+ dt * Bs @ u.iloc[k, :]
# Indoor temperature
y = Cs @ θ + Ds @ u.to_numpy().T
y = y.T
# HVAC heat flow
q_HVAC = Kp * (data['Ti'] - y[:, 0])
# Plot results
return y, q_HVAC, data
dt = 360 # time step, s
y, q_HVAC, data = P_control(
filename, start_date, end_date, dt,
Af, Bf, Cf, Df, Kpf, Tisp)
tuto.plot_results(y, q_HVAC, data)
Figure 1. Simulation in free-running with weather data using Euler explicit method of integration and a time step of 400 s. a) Indoor and outdoor temperatures. b) Solar and HVAC heat flow rates.
Since the building is in free-running, the energy consumption for HVAC is zero, $Q_{HVAC} = 0 \text{ J}$.
The energy consumption over the time range $[0,\, t_{final}]$ is:
$$ Q_{HVAC} = \int_{0}^{t_{final}} \left | q_{HVAC} \right | \,dt$$$$= \Delta t \sum_{k=0}^{n} \left |q_{HVAC(k)} \right| \text{ [J]}$$$$= \frac{1}{3.6·10^{6}} \Delta t \sum_{k=0}^{n}\left |q_{HVAC(k)} \right| \text{ [kW h]}$$where:
- $\Delta t$ is the time step, s.
- $\left |q_{HVAC(k)} \right|$ - absolute values of the vector of heat flow rates, W.
print(f"Q_HVAC = {dt * sum(abs(q_HVAC)) / 1000 / 3600:.1f} kW·h")
Q_HVAC = 0.0 kW·h
Model for building with "perfect" proportional controller¶
A "perfect" proportional controller would have an infinite gain, $K_p \rightarrow \infty$. The change in the controller gain implies a change of the maximum time step.
Kpc = 1e3 # perfect controller Kp -> ∞
TCa = tuto.thermal_circuit(Kpc)
[Ac, Bc, Cc, Dc] = dm4bem.tc2ss(
TCa['A'], TCa['G'], TCa['C'], TCa['b'], TCa['f'], TCa['y'])
dt_max = min(-2. / np.linalg.eig(Ac)[0])
print(f'For Kp = {Kpc:.0e}, the maximum time step is {dt_max:.0f} s')
For Kp = 1e+03, the maximum time step is 59 s
dt = 50 # time step, s
y, q_HVAC, data = P_control(
filename, start_date, end_date, dt,
Ac, Bc, Cc, Dc, Kpc, Tisp)
tuto.plot_results(y, q_HVAC, data)
Figure 2. Simulation of an almost perfect controller with weather data using Euler explicit method of integration and a time step of 50 s. a) Indoor and outdoor temperatures. b) Solar and HVAC heat flow rates.
In the case of an almost perfect controller, the indoor temperature is maintained rather constant and the HVAC heat flow rate $q_{HVAC}$ is positive for heating and negative for cooling.
The energy consumption over the time range $[0, t_{final}]$ is:
$$ Q_{HVAC} = \frac{1}{3.6·10^{6}} \Delta t \sum_{k=0}^{n}\left |q_{HVAC(k)} \right| \text{ [kW h]}$$where:
- $\Delta t$ is the time step, s.
- $q_{(k)HVAC}$ - the values of the vector of heat flow rates, W.
print(f"Q_HVAC = {dt * sum(abs(q_HVAC)) / 1000 / 3600:.1f} kW·h")
Q_HVAC = 168.0 kW·h
Heating and cooling with dead band¶
Maintaing the indoor temperature at a constant level is not energy efficient. It is preferable to have a deadband for the indoor temperature in which the building is in free running and to use the HVAC system when the indoor temperature is out the deadband.
An appraoch for the simulation of heating and cooling with deadband is to use both models, free-running and perfect controlled. In this case, the time step is the minimum between the two time steps.
The code is implemented in function tuto.switch_models() which is very similar to the function tuto.P_control(). The differences consist in:
def switch_models(filename, start_date, end_date, dt,
Af, Bf, Cf, Df, Kpf,
Ac, Bc, Cc, Dc, Kpc,
Tisp, DeltaT):
t, u, data = tuto.inputs(filename, start_date, end_date, dt, Tisp)
# initial values for temperatures
temp_exp = 0 * np.ones([Af.shape[0], u.shape[0]])
Tisp = Tisp * np.ones(u.shape[0])
y = np.zeros(u.shape[0])
y[0] = Tisp[0]
q_HVAC = 0 * np.ones(u.shape[0])
# integration in time
I = np.eye(Af.shape[0])
for k in range(u.shape[0] - 1):
if y[k] < Tisp[k] or y[k] > DeltaT + Tisp[k]:
temp_exp[:, k + 1] = (I + dt * Ac) @ temp_exp[:, k]\
+ dt * Bc @ u.iloc[k, :]
y[k + 1] = Cc @ temp_exp[:, k + 1] + Dc @ u.iloc[k + 1]
q_HVAC[k + 1] = Kpc * (Tisp[k + 1] - y[k + 1])
else:
temp_exp[:, k + 1] = (I + dt * Af) @ temp_exp[:, k]\
+ dt * Bf @ u.iloc[k, :]
y[k + 1] = Cf @ temp_exp[:, k + 1] + Df @ u.iloc[k]
q_HVAC[k + 1] = 0
return y, q_HVAC, data
For the indoor temperature setpoint Tisp and the accepted deadband in temperature DeltaT.
θ_isp = 20 # indoor temperature setpoint, °C
Δθ = 2.5 # deadband of indoor temperature, °C
Large proportional gain for the controller¶
If the controller has a large proportional gain, the temperature is maintained in the limits. However, the heat flow of the HVAC system and the indoor tempertaure vary too quickly, phenomena called pumping (Figure 3). Pumping is different of instability in the sens that the variables have limited values between bounds. Pumping is characterized by repeated oscilation between the bounds of the variables which causes wear and tear of the equipement.
dt = 50
Kpc = 1e3 # no controller Kp -> 0
TCa = tuto.thermal_circuit(Kpc)
[Ac, Bc, Cc, Dc] = dm4bem.tc2ss(
TCa['A'], TCa['G'], TCa['C'], TCa['b'], TCa['f'], TCa['y'])
y, q_HVAC, data = switch_models(
filename, start_date, end_date, dt,
Af, Bf, Cf, Df, Kpf,
Ac, Bc, Cc, Dc, Kpc,
Tisp, Δθ)
tuto.plot_results(y, q_HVAC, data)
ax = plt.gca()
ax.set_ylim([0, 2000])
(0.0, 2000.0)
Figure 3. Simulation of heating and cooling with deadband by using a high gain controller. a) Indoor and outdoor temperature. b) Solar and HVAC heat flow rates.
print(f"Q_HVAC = {dt * sum(abs(q_HVAC)) / 1000 / 3600:.1f} kW·h")
Q_HVAC = 179.9 kW·h
Large proportional gain for the controller with reduced time step¶
dt = 25
Kpc = 1e3 # no controller Kp -> 0
TCa = tuto.thermal_circuit(Kpc)
[Ac, Bc, Cc, Dc] = dm4bem.tc2ss(
TCa['A'], TCa['G'], TCa['C'], TCa['b'], TCa['f'], TCa['y'])
y, q_HVAC, data = switch_models(
filename, start_date, end_date, dt,
Af, Bf, Cf, Df, Kpf,
Ac, Bc, Cc, Dc, Kpc,
Tisp, Δθ)
tuto.plot_results(y, q_HVAC, data)
ax = plt.gca()
ax.set_ylim([0, 1000]);
Figure 4. Heating and cooling with deadband, high Kp and small time step. a) Indoor and outdoor temperature. b) Solar and HVAC heat flow rates.
Note that the energy consumption depends on the time step.
print(f"Q_HVAC = {dt * sum(abs(q_HVAC)) / 1000 / 3600:.1f} kW·h")
Q_HVAC = 114.6 kW·h
Optimal controller gain¶
The instability or pumping effect may be reduced by using a smaller gain for the controller.
dt = 50 # time step, s
Kpc = 1.5e2 # controller Kp
TCa = tuto.thermal_circuit(Kpc)
[Ac, Bc, Cc, Dc] = dm4bem.tc2ss(
TCa['A'], TCa['G'], TCa['C'], TCa['b'], TCa['f'], TCa['y'])
y, q_HVAC, data = switch_models(
filename, start_date, end_date, dt,
Af, Bf, Cf, Df, Kpf,
Ac, Bc, Cc, Dc, Kpc,
Tisp, Δθ)
tuto.plot_results(y, q_HVAC, data)
ax = plt.gca()
ax.set_ylim([0, 1000]);
Figure 5. Heating and cooling with dead-band and a sub-optimal gain for the controller. a) Indoor and outdoor temperature. b) Solar and HVAC heat flow rates.
The energy consumption is significantly smaller as compared with the "perfect" controller.
print(f"Q_HVAC = {sum(dt * abs(q_HVAC)) / 1000 / 3600:.1f} kW·h")
Q_HVAC = 119.5 kW·h
Only heating¶
If there is only heating, then the HVAC system controls the indoor temperature when it is smaller than the setpoint and leaves the building in free running if the indor temperature is larger than the setpoint.
Heating is implemented in tuto.heat() function, which is similar to switch_models() fuction with the exception of the numerical integration.
t, u, data = tuto.inputs(filename, start_date, end_date, dt, Tisp)
# initial values for temperatures
temp_exp = 0 * np.ones([Af.shape[0], u.shape[0]])
Tisp = Tisp * np.ones(u.shape[0])
y = np.zeros(u.shape[0])
y[0] = Tisp[0]
q_HVAC = 0 * np.ones(u.shape[0])
I = np.eye(Af.shape[0])
for k in range(u.shape[0] - 1):
if y[k] < data['Ti'][k]:
temp_exp[:, k + 1] = (I + dt * Ac) @ temp_exp[:, k]\
+ dt * Bc @ u.iloc[k, :]
y[k + 1] = Cc @ temp_exp[:, k + 1] + Dc @ u.iloc[k]
q_HVAC[k + 1] = Kpc * (data['Ti'][k + 1] - y[k + 1])
else:
temp_exp[:, k + 1] = (I + dt * Af) @ temp_exp[:, k]\
+ dt * Bf @ u.iloc[k, :]
y[k + 1] = Cf @ temp_exp[:, k + 1] + Df @ u.iloc[k]
q_HVAC[k + 1] = 0
tuto.plot_results(y, q_HVAC, data)
Figure 6. Heating. a) Indoor and outdoor temperature. b) Solar and HVAC heat flow rates.
If only heating is used, the energy consumption is smaller.
print(f"Q_HVAC = {dt * sum(abs(q_HVAC)) / 1000 / 3600:.1f} kW·h")
Q_HVAC = 69.1 kW·h
However, there is overheating. This could be aleviate by controlling the solar radiation and the free-cooling by ventilation modelled as inputs.