A Review of Mathematical Models of Building Physics and Energy Technologies for Environmentally Friendly Integrated Energy Management Systems

Abstract

The Energy Management System (EMS) is an efficient technique to monitor, control and enhance the building performance. In the state-of-the-art, building performance analysis is separated into building simulation and control management: this may cause inaccuracies and extra operating time. Thus, a coherent framework to integrate building physics with various energy technologies and energy control management methods is highly required. This framework should be formed by simplified but accurate models of building physics and building energy technologies, and should allow for the selection of proper control strategies according to the control objectives and scenarios. Therefore, this paper reviews the fundamental mathematical modeling and control strategies to create such a framework. The mathematical models of (i) building physics and (ii) popular building energy technologies (renewable energy systems, common heating and cooling energy systems and energy distribution systems) are first presented. Then, it is shown how the collected mathematical models can be linked. Merging with two frequently used EMS strategies, namely rule-based and model predictive controls, is discussed. This work provides an extendable map to model and control buildings and intends to be a foundation for building researchers, designers and engineers.

1. Introduction

The building sectors account for around 36% of energy consumption and 39% of carbon emissions of the globe [1]. Thus, improving the building energy performance by increasing the use of clean and renewable energy systems is a promising direction to realize resource conservation and face climate change. With the development of the concept of ’nearly zero-energy building’, some novel technologies, such as renewable energy resources, onsite storage systems and environmentally friendly devices in buildings, have gained more interest [2].

As an alternative to producing electricity from fossil fuels, renewable energy resources based on solar and wind provide clean and low-cost energy for buildings and have been increasingly used in recent years (Photovoltaic Panels (PVs), Solar Thermal Panels (STPs), Wind Turbines (WTs), etc.) [3]. The existing literature shows the generation model of solar and wind power under different weather conditions [4,5], as well as algorithms to maximize their usage ratio in buildings [6,7].

However, the inherent uncertainty and intermittence of renewable generation bring new challenges to the power systems. This means that the onsite renewable generation profile may not be consistent with the consumption profile in buildings. Hence, some energy storage technologies (electricity storage or thermal storage) are often used simultaneously with renewable resources as effective techniques to deal with the mismatch between renewable generation and local consumption [8]. Meanwhile, energy storage technology can improve energy flexibility and decrease the operation cost [9]. It is shown that a combination of renewable generation and energy storage systems can reduce up to 30% of the annual cost and increase up to 29% of the self-consumption proportion of renewable energy [10].

Environmentally friendly technologies with a high Coefficient of Performance (COP)—for example, heat pumps and Combined Heating and Power (CHP)—have also been increasingly studied and developed [11,12]. These devices can be combined with water-based energy distribution systems that allow for low-temperature heating [13] and high-temperature cooling [14] to cover the heating and cooling demand of the buildings.

With so many technologies co-existing in one building, a smart Energy Management System (EMS) solution to operate and make them profitable for users needs to be considered. There are generally two kinds of EMS: passive EMS and active EMS [15]. The former improves the energy efficiency in buildings by predicting the load curves and influencing users’ awareness [16], whereas the latter monitors and controls the elements of building energy systems via a set of well-designed control strategies [17]. EMS has been researched and applied to span from the device level [18] and the single building level [19] up to the building community level [20]. Typical EMS implementations for buildings include defining a series of control rules dependent on the energy or electricity price thresholds [21], executing mixed integer optimization under time-of-use electricity price and incentive policy [22], establishing a game theory model between user behaviours and the electricity price policy [23] and using a bionic algorithm, such as genetic algorithms [24] or particle swarm optimization [19], to optimize the operation of household appliances. The EMS can also interact with the grid side to determine the incentives [25] or reduce the energy bill [26]. For the grid operators, this interaction helps to protect the power grid from the risk of power blackouts by reducing the peak-to-average ratio [27].

1.1. Research Motivation: Integration

Studying complex EMSs in buildings requires the use of building simulation and programming tools separately or collectively. Building simulation tools are mature techniques to analyze the building performance under different targets; for example, energy consumption and thermal comfort. Several building simulation tools (e.g., TRNSYS, EnergyPlus, IDA-ICE) are widely developed to provide a precise model for the buildings and technologies [11,28,29,30,31]. However, all building simulation tools face a similar situation when integrated with an EMS: they are only able to account for the modeling of buildings (and sometimes also a few standard equipment), but they cannot be merged in a seamless way with the EMS based on different control strategies. As a matter of fact, usually the data from simulation tools must be collected into another programming tool (e.g., MATLAB, Python) to complete the control part. This difficulty becomes more apparent when the prediction is taken into account in an EMS. For example, an EMS covering energy flexibility conversion, routing and storage options was investigated in [32], where the modeling and control parts were carried out in TRNSYS and MATLAB separately. In addition, different EMSs for a community of buildings were studied in [11], sorted out into two parts: first, calculating each building heating energy demand by building simulation tools (TRNSYS or IDA ICE); then, using the data as input values into the programmed EMS in MATLAB. This strategy to undertake the simulation tools for building modeling and the programming tools for executing EMS strategies separately brings difficulties and challenges in setting up a framework, which increases the run time and may cause inaccuracies through the final results. This may especially be inconvenient when exploring an analysis of complex scenarios (such as buildings integrated with energy technologies and EMS).

1.2. Related Works

Different models of energy technologies have been proposed in the literature, cf.

Table 1. Studied technologies in buildings.

Technologies	Models	Applications in Buildings	Main Findings
Photovoltaic Panel (PV)	Analytical [4] Empirical [35]	[6,36,37]	PVs can produce renewable electricity and reduce over 40% of the energy demand due to their passive benefits.
Solar Thermal Panel (STP)	Analytical [38] Empirical [39]	[40,41]	STPs provide renewable heat generation, with overall efficiency ranges from 24–28% (single pass) or 32–34% (double pass).
Wind Turbine (WT)	Analytical [5,42] Empirical [43]	[26]	WTs can produce renewable electricity and are promising in areas with strong wind currents. It grows rapidly worldwide, with an over 11% growth rate.
Combined (Cooling), Heating and Power (C(C)HP)	Analytical [44,45,46] Empirical [47,48]	[41,49,50]	CHP and CCHP can increase the electricity and heat production efficiency from 60% in traditional ways to 90%.
Battery	Analytical [36] Empirical [51]	[20,52,53]	Batteries can deal with the intermittence of renewable generation and shave peak demand, contributing to balancing the demand and supply and maintaining the grid stability, and can achieve an over 20% cost reduction.
Hot Water Storage Tank (HWST)	Analytical [54,55] Empirical [56]	[11,57,58]	HWSTs can shift thermal loads to off-peak hours and increase energy flexibility, providing up to a 7.5% energy saving potential and 5.5% emission reduction potential.
Phase Change Material (PCM)	Analytical [59,60] Empirical [61]	[62,63,64]	PCMs have large thermal storage abilities due to their high latent heat, leading to the same energy saving with less materials.
Air/Ground Source Heat Pump (ASHP/GSHP)	Analytical [65,66] Empirical [67,68,69]	[12,70,71,72]	Heat pumps have a high operating efficiency and low operating cost compared with conditional heating devices, with a COP of up to 5.
Direct/Indirect Evaporative Cooling (DEC/IEC)	Analytical [73,74] Empirical [75]	[76,77]	DEC and IEC are efficient in hot and humid climates, and can decrease the air temperature to wet bulb temperature or dew point temperature.
Absorption/ Compression Chiller (AC/CC)	Analytical [78] Empirical [79,80]	[81,82,83]	The COPs of chillers range from 0.6 to 6, and the consumed energy can be supplied by renewable heat/electricity production. The COP of CCs is usually higher than ACs.

. However, coherent framework that can integrate all of these models in a unified way is still missing. The aforementioned framework should describe how the different technologies interact with each other, and how they can be operated by the EMS (e.g., by selecting certain thresholds or a certain power). A possible solution to the issue of integration is to merge physics-based mathematical models to create comprehensive building energy systems [33]. Then, corresponding EMSs can be easily coded and linked to such building energy systems models. Towards this purpose, bottom-up load models of the residential space heating and cooling load, Domestic Hot Water (DHW) load, appliances and equipment and electric vehicles have been developed and validated in [33]. However, these models were based on traditional heating technologies; as a result, this approach does not cover renewable energy resources and efficient environmentally friendly technologies, such as heat pumps. Another EMS framework focused on microgrid control and optimization based on objectives, constraints and optimization methods was reviewed in [34], but most popular and environmentally friendly technologies and energy distribution systems were not discussed.

Although the techniques for modeling HVAC systems were reviewed in [84,85], it is unclear how the component models can be connected to model the whole energy systems in buildings. In fact, some studies have shown that the integration between modeling and control requires dedicated mathematical models, e.g., modular approaches to connect different devices [86,87]. Another example is the integrated HVAC and DHW production systems discussed in [88], which are given without presenting the models. It is mentioned that there is a great challenge to minimize the mismatching between the design and real energy performance, which highlighted the requirement for unified modeling methods for multi-energy systems. Taking into account the aforementioned concerns, setting up a whole framework that includes both the mathematical models and control strategies for buildings installed with environmentally friendly technologies is still an open problem. Most importantly, this integration should fit different configurations in different buildings.

1.3. Contribution of This Study

This study contributes computationally efficient models suitable for the EMS design for buildings. The aim is to describe a selective and integrated approach. The mathematical models presented in this work can be connected to construct different energy systems in different buildings. In this way, the models of building physics and energy systems cover all types of technologies popularly installed in buildings. The experiments to verify the selected models are not presented here, but they are referred to the validations in the reviewed studies; for example, by comparing with the results of building simulation tools [56] or measurements data [74,89,90]. The discussed models are expandable and integrable. Therefore, a direction is to combine them and to analyze the building level performance. Besides, this paper provides a whole map to establish EMS in buildings based on the different layout and installed technologies, and contributes to the further research of EMS design for researchers, designers and engineers.

The rest of this paper is organized as follows. Section 2, Section 3 and Section 4 present a review of the assemblable mathematical models for building physics, heating and cooling energy systems and onsite energy generation and storage. Section 5 discusses two common EMS strategies, namely Rule-Based Control (RBC) and Model Predictive Control (MPC). Section 6 summarizes the framework and gives concluding remarks.

2. Building Physics

Models for building physics can be categorized as white-box, black-box and grey-box [91]. The first approach is adopted by most building simulation tools, such as TRNSYS and EnergyPlus: this approach consists of estimating the demand and consumption of a specific building based on weather data, building envelope, installed technologies and consumption profiles totally based on the physical relationships. Examples of black-box models are neural networks or autoregressive moving average models, in which, the mathematical relationships between inputs (e.g., weather data) and outputs (e.g., energy consumption) are purely based on collected data. A grey-box method is a hybrid approach combining white-box and black-box models, since it is based on physical formulations (e.g., heat transfer laws) with parameters tuned by measured data.

A detailed review of white-box modeling and black-box modeling was presented in [84,85]. These reviews mention that the white-box models are often developed based on a number of assumptions and parameters, which causes a complex construction process and low accuracy. On the other hand, the accuracy of black-box models usually depends on the types of the selected methods. A large number of training data are required to develop the black-box models, with a relatively lower generalization capability. This means that the generated relationships are valid for a specific building but unsuitable for other different building configurations. To achieve a framework that takes into account both building physics and heterogeneous energy technologies, the grey-box approach based on the widely used thermal Resistance and Capacitance (RC) network model is one of the most adopted approaches. An RC model exploits the analogy between the energy flow and electric flow, since a capacitor represents an accumulation of energy, whereas the resistance represents a dissipation of energy. The two main factors that influence the accuracy and complexity of an RC network are its spatial resolution and order. Spatial resolution simply refers to the option of aggregating all rooms as a single zone, or describing each room with a separate RC model. The order of an RC model refers to the number of lumped elements in the model; for example, one capacitor (1C) [92] or two capacitors (2C) [93,94]. The second option is more capable of describing phenomena such as the thermal capacity of the air and building fabric [95]. According to their order, the most standard RC models are the so-called 2R-1C [96] or 3R-2C [97]. The 3R-2C model is usually suggested as a trade-off between complexity and accuracy, and examples reported in the literature include a large retail store building with 28 zones in California [98] and an office building with 7 zones and 33 surfaces [99], indicating that this model can be adjusted for buildings with a wide range of sizes.

Apart from the energy flow and heat exchange, infiltration due to leaks and cracks in the building envelope has a non-negligible influence on the total energy demand [100]. Measuring the infiltration rate in buildings is essential to control and manage the energy consumption and improve indoor thermal levels. Various methods were developed to achieve this value [101]. For simplicity, the infiltration can be assessed by the airtightness, i.e., the leakage–infiltration ratio of buildings [102]. This value expresses a linear relationship between the air leakage rate and the infiltration rate, and is considered over a period of time.

The schematic representation of a general multi-zone 3R-2C network is shown in Figure 1, with each node in the model representing a point of interest. Once the layout of one floor in the building is determined, the multi-zone 3R-2C network is expressed as follows:

$$\frac{1}{2}{M}_{\mathrm{exwall},i}{C}_{\mathrm{exwall}}{\dot{T}}_{\mathrm{exwall},i}^{\mathrm{in}}=h_{\mathrm{in}}{A}_{\mathrm{exwall},i}(T_{\mathrm{zone},i}-T_{\mathrm{exwall},i}^{\mathrm{in}})+\frac{T_{\mathrm{exwall},i}^{\mathrm{out}}-T_{\mathrm{exwall},i}^{\mathrm{in}}}{R_{\mathrm{exwall}}},$$

(1)

$$\frac{1}{2}{M}_{\mathrm{exwall},i}{C}_{\mathrm{exwall}}{\dot{T}}_{\mathrm{exwall},i}^{\mathrm{out}}=h_{\mathrm{out}}{A}_{\mathrm{exwall},i}(T_{\mathrm{amb}}-T_{\mathrm{exwall},i}^{\mathrm{out}})+\frac{T_{\mathrm{exwall},i}^{\mathrm{in}}-T_{\mathrm{exwall},i}^{\mathrm{out}}}{R_{\mathrm{exwall}}},$$

(2)

$$\begin{array}{cc}\hfill M_{\mathrm{inwall},i-j}{C}_{\mathrm{inwall}}{\dot{T}}_{\mathrm{inwall},i-j}=& h_{\mathrm{in}}{A}_{\mathrm{inwall},i-j}(T_{\mathrm{zone},i}-T_{\mathrm{inwall},i-j})\hfill \\ \hfill +& h_{\mathrm{in}}{A}_{\mathrm{inwall},i-j}(T_{\mathrm{zone},j}-T_{\mathrm{inwall},i-j}),\hfill \end{array}$$

(3)

$$\begin{array}{cc}\hfill M_{\mathrm{air},i}{C}_{\mathrm{air}}{\dot{T}}_{\mathrm{zone},i}=& h_{\mathrm{in}}{A}_{\mathrm{exwall},i}(T_{\mathrm{exwall},i}^{\mathrm{in}}-T_{\mathrm{zone},i})+\frac{T_{\mathrm{amb}}-T_{\mathrm{zone},i}}{R_{\mathrm{win}}}+{\dot{m}}_{\inf }{C}_{\mathrm{air}}(T_{\mathrm{amb}}-T_{\mathrm{zone},i})\hfill \\ \hfill +& \sum _j{h}_{\mathrm{in}}{A}_{\mathrm{inwall},i-j}(T_{\mathrm{inwall},i-j}-T_{\mathrm{zone},i})+Q_{\mathrm{HVAC},i}+Q_{\mathrm{int},i}+Q_{\mathrm{solar},i}.\hfill \end{array}$$

(4)

This model is applied for a variety of buildings and the coefficients can be calculated based on the information of the building envelope. It can also be reduced to a single-zone model by assuming $n=1$.

3. Onsite Energy Generation and Storage

This section concludes a review of the mathematical models of onsite energy generation (PV, STP, WT, CHP) and energy storage units (battery, HWST, PCM).

3.1. Photovoltaic Panel

A PV produces electricity based on solar energy. The generated power is influenced by the weather conditions, among which, the ambient temperature and solar radiation are the main factors. For simplicity, most models assume that the total radiation is falling on the solar array, and the angle of incidence is not considered [26]. A linear model used by [37] expresses the generated power $P_{\mathrm{PV}}$ (W) as:

$$P_{\mathrm{PV}}={\eta }_{\mathrm{PV}}{A}_{\mathrm{PV}}G\left(1-0.005(T_{\mathrm{amb}}-25)\right),$$

(5)

where ${\eta }_{\mathrm{PV}}$ (%) indicates the conversion efficiency, which is typically equal to 16%, and G (W/m${}^2$) is the solar radiation. An alternative model to (5) is proposed in [6]:

$$P_{\mathrm{PV}}=\frac{G}{1000}\left(1+{\alpha }_{\mathrm{cur}}\Delta T\right)\left(1+{\alpha }_{\mathrm{vol}}\Delta T\right)P_{\mathrm{PV},\mathrm{rated}},$$

(6)

$$\Delta T=T_{\mathrm{cell}}-25,$$

(7)

where ${\alpha }_{\mathrm{cur}}$ ($\%/{}^{\circ }\mathrm{C}$) and ${\alpha }_{\mathrm{vol}}$ ($\%/{}^{\circ }\mathrm{C}$) are the temperature sensitivity of the photovoltaic generated current and voltage, their typical values are ${\alpha }_{\mathrm{cur}}=0.05$%/${}^{\circ }\mathrm{C}$ and ${\alpha }_{\mathrm{vol}}=-0.33$$\%/{}^{\circ }\mathrm{C}$, $P_{\mathrm{PV},\mathrm{rated}}$ (W) is the nominal power of the PV in the standard test conditions ($G=1000$ W/m${}^2$, $T_{\mathrm{cell}}=25{}^{\circ }\mathrm{C}$, AM = 1.5) and $T_{\mathrm{cell}}$ (${}^{\circ }\mathrm{C}$) is the cell temperature, which can be estimated by [41]:

$$T_{\mathrm{cell}}=T_{\mathrm{amb}}+\frac{1}{32}G.$$

(8)

The model (5) is a linear model, whereas (6)–(8) constitute a quadratic one. It is possible to tune the parameters of both models in such a way that their results are similar in a certain operating region. However, the linear model is able to estimate the generation for a wider range of cell temperatures, whereas the quadratic one predicts the generation around the nominal temperature more accurately. A combination of them can be used, e.g., the quadratic model is used when $T_{\mathrm{amb}}$ is the nominal temperature and the linear model is used when $|\Delta T|$ is large.

3.2. Solar Thermal Panel

A STP produces thermal heat with solar energy. Solar radiation and ambient temperature are also the dominant factors of the STP performance. The usable heat $Q_{\mathrm{STP}}$ (W) produced by a flat plate STP can be calculated with the Hottel–Whillier–Bliss equation [103]:

$$Q_{\mathrm{STP}}=A_{\mathrm{STP}}{F}_{\mathrm{R}}\left[\tau G-U_{\mathrm{STP},\mathrm{L}}(T_{\mathrm{STP},\mathrm{in}}-T_{\mathrm{amb}})\right]={\dot{m}}_{\mathrm{STP}}{C}_{\mathrm{w}}(T_{\mathrm{STP},\mathrm{out}}-T_{\mathrm{STP},\mathrm{in}}),$$

(9)

where $A_{\mathrm{STP}}$ (m${}^2$) refers to the area of the collector, $F_{\mathrm{R}}$ (%) is the heat removal efficiency factor, $\tau$ (%) is the transmittance absorptivity, $U_{\mathrm{STP},\mathrm{L}}$ (W/(m${}^2$ K)) is the heat loss coefficient, ${\dot{m}}_{\mathrm{STP}}$ (kg/s) is the water mass flow rate, $T_{\mathrm{STP},\mathrm{in}}$ (${}^{\circ }\mathrm{C}$) and $T_{\mathrm{STP},\mathrm{out}}$ (${}^{\circ }\mathrm{C}$) is the inlet and outlet water temperature and $C_{\mathrm{w}}$ (J/(kg K)) is the specific heat of water. It is assumed that $F_{\mathrm{R}}$, $\tau$, $U_{\mathrm{STP},\mathrm{L}}$ are constants for a given collector, water flow rate (typically 0.005 kg/(s m${}^2$)) and fluid inlet temperature [104]. In particular, the technical characteristics are assumed to be $F_R=0.67$%, $\tau =0.77$%, $U_{\mathrm{STP},\mathrm{L}}=7.47$ W/(m${}^2$ K), which are typical values for double-glazed flat-plate collectors [105]. The model (9) calculates the usable heat produced by STP based on a function of its area and solar radiation, and is simplified by assuming the coefficients as constants.

The combination of PV and STP forms a new technology called photovoltaic thermal collectors, which simultaneously converts solar energy into heat and electricity [106]. It is noted that the electricity and heat production efficiency of photovoltaic thermal collectors is lower than that of the separate PV or STP, but that its overall efficiency is higher [38].

3.3. Wind Turbine

A WT is a local generation device based on wind energy. It works when the wind speed $v_{\mathrm{wind}}$ (m/s) is within a certain range (larger than cut-in velocity $v_{\mathrm{c}}$ (m/s) and less than cut-off velocity $v_{\mathrm{f}}$ (m/s)). In fact, the power $P_{\mathrm{WT}}$ (W) generated by the WT can be calculated by the following piece-wise function of the wind velocity [107]:

$$P_{\mathrm{WT}}=\left\{\begin{array}{cc}0\hfill & \mathrm{if}{v}_{\mathrm{f}}\le v_{\mathrm{wind}}\mathrm{or}{v}_{\mathrm{wind}}\le v_{\mathrm{c}}\hfill \\ \frac{v_{\mathrm{wind}}^3-v_c^3}{v_{\mathrm{rated}}^3-v_{\mathrm{c}}^3}{P}_{\mathrm{WT},\mathrm{rated}}\hfill & \mathrm{if}{v}_{\mathrm{c}}<v_{\mathrm{wind}}<v_{\mathrm{rated}}\hfill \\ P_{\mathrm{WT},\mathrm{rated}}\hfill & \mathrm{if}{v}_{\mathrm{rated}}\le v_{\mathrm{wind}}<v_f\hfill \end{array}\right.,$$

(10)

where $P_{\mathrm{WT},\mathrm{rated}}$ (W) is the rated electrical power and $v_{\mathrm{rated}}$ (m/s) is the rated wind speed. Other approaches calculate $P_{\mathrm{WT}}$ in an alternative but analogous way [26]:

$$P_{\mathrm{WT}}=\left\{\begin{array}{cc}0\hfill & \mathrm{if}{v}_{\mathrm{f}}\le v_{\mathrm{wind}}\mathrm{or}{v}_{\mathrm{wind}}\le v_{\mathrm{c}}\hfill \\ 1/2{\rho }_{\mathrm{air}}\pi R_{\mathrm{WT}}^2{v}_{\mathrm{wind}}^3{C}_{\mathrm{WT}}\hfill & \mathrm{if}{v}_{\mathrm{c}}<v_{\mathrm{wind}}<v_{\mathrm{rated}}\hfill \\ P_{\mathrm{WT},\mathrm{rated}}\hfill & \mathrm{if}{v}_{\mathrm{rated}}\le v_{\mathrm{wind}}<v_{\mathrm{f}}\hfill \end{array}\right.,$$

(11)

where ${\rho }_{\mathrm{air}}$ ($\mathrm{kg}/{\mathrm{m}}^3$) is the air density, $R_{\mathrm{WT}}$ (m) is the blade radius and $C_{\mathrm{WT}}$ (%) is the power coefficient. It is assumed that ${\rho }_{\mathrm{air}}=1.1839\mathrm{kg}/{\mathrm{m}}^3$, which is the air density at sea level and 25 ${}^{\circ }\mathrm{C}$ and $C_{\mathrm{WT}}=0.4$%, which is a typical value in the literature. These two formulations represent a simple way to integrate data of wind velocity in the generation of wind energy.

3.4. Combined Heating and Power

CHP or micro-CHP is a class of equipment in which the waste heat of an electrical generator is collected to cover thermal loads [47]. A CHP system usually consists of the following parts: (1) prime mover, e.g., combustion turbine, engine or microturbine; (2) electric generator; (3) heat recovery unit, which transforms wasted heat to reusable thermal energy.

An analytical dynamic model for a micro-CHP with a Stirling engine taking into consideration the start-up and cool-down periods and partial load performance was developed in [108], and an experimentally validated model of an internal-combustion-engines-driven, natural-gas-fueled device based on a product of Marathon engine systems ecopower was studied and tested in [47], demonstrating that the generated heat and electricity power were nearly linear functions of the engine speed within the operating range. Another empirical model for the engine was polynomials presented in [41], and is suitable for a wider range of CHP rated power. According to this model, the electricity and thermal heat generation efficiency is estimated by:

$${\epsilon }_{\mathrm{CHP},\mathrm{e}}=0.4679{\eta }_{\mathrm{CHP}}^3-1.1705{\eta }_{\mathrm{CHP}}^2+0.998{\eta }_{\mathrm{CHP}}+0.002,$$

(12)

$${\epsilon }_{\mathrm{CHP},\mathrm{h}}=-0.4845{\eta }_{\mathrm{CHP}}^3+1.0934{\eta }_{\mathrm{CHP}}^2-0.8379{\eta }_{\mathrm{CHP}}+0.7221,$$

(13)

where ${\eta }_{\mathrm{CHP}}$ (%) is the part load ratio calculated by:

$${\eta }_{\mathrm{CHP}}=P_{\mathrm{CHP}}/P_{\mathrm{CHP},\mathrm{rated}},$$

(14)

where $P_{\mathrm{CHP}}$ (W) and $P_{\mathrm{CHP},\mathrm{rated}}$ (W) are the generated and rated electricity power. Then, the relationships between the fuel input $Q_{\mathrm{fuel}}$ (W) and electricity and thermal outputs are:

$$P_{\mathrm{CHP}}={\epsilon }_{\mathrm{CHP},\mathrm{e}}{Q}_{\mathrm{fuel}},$$

(15)

$$Q_{\mathrm{CHP}}={\epsilon }_{\mathrm{CHP},\mathrm{h}}{Q}_{\mathrm{fuel}}.$$

(16)

The model directly focuses on the relationship between fuel consumption and produced energy, The overall system efficiency is usually used to assess the performance of a CHP, which is the sum of produced electricity and thermal heat divided by the total fuel energy input [109].

3.5. Battery

A battery provides the potential to purchase or store electricity during off-peak hours and to sell or use electricity during peak hours in buildings. The models for batteries derived from the literature can be classified into linear and nonlinear ones [110]. The recommended model in this framework is the linear one [53]:

$$0\le {\dot{Q}}_{\mathrm{ch}}\le {\dot{Q}}_{\mathrm{ch}}^{\max },0\le {\dot{Q}}_{\mathrm{dis}}\le {\dot{Q}}_{\mathrm{dis}}^{\max },$$

(17)

$${\dot{Q}}_{\mathrm{ch}}{\dot{Q}}_{\mathrm{dis}}=0,$$

(18)

$$\dot{S}=-{\eta }_{\mathrm{loss}}S+{\eta }_{\mathrm{ch}}{\dot{Q}}_{\mathrm{ch}}-{\eta }_{\mathrm{dis}}{\dot{Q}}_{\mathrm{dis}},$$

(19)

$$S_{\min }\le S\le S_{\max },$$

(20)

where ${\dot{Q}}_{\mathrm{ch}}$ (W/s) and ${\dot{Q}}_{\mathrm{dis}}$ (W/s) are the charging and discharging rates, and they cannot be positive simultaneously for the safe operation and long life span of the battery. ${\dot{Q}}_{\mathrm{ch}}^{\max }$ (W/s) and ${\dot{Q}}_{\mathrm{dis}}^{\max }$ (W/s) are the limited charging and discharging rates for the safe operation of the battery, S (W) is the state of charge, ${\eta }_{\mathrm{loss}}$ (%) is the loss proportion due to energy leakage and ${\eta }_{\mathrm{ch}}$ (%) and ${\eta }_{\mathrm{dis}}$ (%) denote the charging and discharging efficiency. Typically, ${\eta }_{\mathrm{loss}}\ll 1$, ${\eta }_{\mathrm{ch}}<1<{\eta }_{\mathrm{dis}}$ and $S_{\min }$ and $S_{\max }$ are usually 20% and 80% of the total battery capacity. This model is easy to implement compared with nonlinear models: it only requires knowledge of the initial state of charge of the battery [34].

3.6. Hot Water Storage Tank

The HWST is a cost-effective technology for thermal heat storage. Tanks based on stratification give a higher efficiency than a fully mixed tank of the same size [111]. In this work, a stratified HWST that enables a high temperature layer at the top and a low temperature layer at the bottom of the tank is used.

The modeling methods of a stratified HWST can be categorized in terms of its dimension. Two-dimensional [112,113] or three-dimensional [114,115] solutions based on Computational Fluid Dynamics (CFD) techniques are accurate but require long calculations. As a result, the one-dimensional tank models are highly regarded, and some studies have shown that their accuracy can be comparable to CFD models [54]. Different approaches and configurations have been proposed for one-dimensional modelings, such as incorporating inflow or outflow buoyancy and a mixing process [55], a node-mixing model for a pressurized water tank with two immersed heat exchangers [116], a mixing factor parameter determined by experimental measurements [56] or including PCM to improve the storage capacity of the tank [89].

The schematic representation of a stratified HWST with five nodes (one node represents one layer) is presented in Figure 2. Different energy technologies can be integrated with it. The energy balance equation in each layer is expressed as follows:

$$\left[\begin{array}{cc}\hfill & \mathrm{Evolution}\mathrm{of}\hfill \\ \hfill & \mathrm{thermal}\hfill \\ \hfill & \mathrm{energy}\hfill \end{array}\right]=\left[\begin{array}{cc}\hfill & \mathrm{Thermal}\mathrm{heat}\hfill \\ \hfill & \mathrm{provided}\mathrm{by}\hfill \\ \hfill & \mathrm{devices}\hfill \end{array}\right]+\left[\begin{array}{cc}\hfill & \mathrm{Heat}\mathrm{transfer}\hfill \\ \hfill & \mathrm{through}\hfill \\ \hfill & \mathrm{water}\mathrm{flow}\hfill \end{array}\right]+\left[\begin{array}{cc}\hfill & \mathrm{Vertical}\mathrm{heat}\hfill \\ \hfill & \mathrm{conduction}\hfill \end{array}\right]+\left[\begin{array}{cc}\hfill & \mathrm{Heat}\mathrm{transfer}\hfill \\ \hfill & \mathrm{with}\mathrm{tank}\hfill \\ \hfill & \mathrm{surface}\hfill \end{array}\right].$$

The size and configuration of this tank is similar to that presented in [11]. The top layer provides DHW for the occupants, and the lower layer is connected to Heating, Ventilation and Air Conditioning (HVAC) for space heating. Technologies such as GSHP, ASHP and micro-CHP are installed in nodes 1–4. An electric heater is used in node 5 to ensure the minimum temperature of DHW, which is typically 55 ${}^{\circ }\mathrm{C}$; for example, in the Finnish building regulations [117]. Domestic Cold Water (DCW) or city water that is around 15 ${}^{\circ }\mathrm{C}$ enters the bottom layer with the same speed of DHW. Let $T_i$ (${}^{\circ }\mathrm{C}$), $i=1,2\cdots 5$ denote the temperatures of five layers from the bottom to the top. The analytical model based on energy balance equations validated by [55,89] for HWST is written as follows:

$$\begin{array}{cc}\hfill M_1{C}_{\mathrm{w}}{\dot{T}}_1=& S_{\mathrm{CHP}}{Q}_{\mathrm{CHP}}/4+{\dot{m}}_{\mathrm{DHW}}{C}_{\mathrm{w}}{T}_{DCW}+{\dot{m}}_{\mathrm{SH}}{C}_{\mathrm{w}}{T}_{\mathrm{SH}}+({\dot{m}}_{\mathrm{GSHP}}+{\dot{m}}_{\mathrm{ASHP}})C_{\mathrm{w}}{T}_2\hfill \\ \hfill & -({\dot{m}}_{\mathrm{DHW}}+{\dot{m}}_{\mathrm{SH}}+{\dot{m}}_{\mathrm{GSHP}}+{\dot{m}}_{\mathrm{ASHP}})C_{\mathrm{w}}{T}_1+Q_1+Q_{\mathrm{case},1},\hfill \end{array}$$

(21)

$$\begin{array}{cc}\hfill M_2{C}_{\mathrm{w}}{\dot{T}}_2=& S_{\mathrm{CHP}}{Q}_{\mathrm{CHP}}/4+({\dot{m}}_{\mathrm{DHW}}+{\dot{m}}_{\mathrm{SH}})C_{\mathrm{w}}{T}_1+({\dot{m}}_{\mathrm{GSHP}}+{\dot{m}}_{\mathrm{ASHP}})C_{\mathrm{w}}{T}_3\hfill \\ \hfill & -({\dot{m}}_{\mathrm{DHW}}+{\dot{m}}_{\mathrm{SH}}+{\dot{m}}_{\mathrm{GSHP}}+{\dot{m}}_{\mathrm{ASHP}})C_{\mathrm{w}}{T}_2+Q_2+Q_{\mathrm{case},2},\hfill \end{array}$$

(22)

$$\begin{array}{cc}\hfill M_3{C}_{\mathrm{w}}{\dot{T}}_3=& S_{\mathrm{CHP}}{Q}_{\mathrm{CHP}}/4+({\dot{m}}_{\mathrm{DHW}}+{\dot{m}}_{\mathrm{SH}})C_{\mathrm{w}}{T}_2+({\dot{m}}_{\mathrm{GSHP}}+{\dot{m}}_{\mathrm{ASHP}})C_{\mathrm{w}}{T}_4\hfill \\ \hfill & -({\dot{m}}_{\mathrm{DHW}}+{\dot{m}}_{\mathrm{SH}}+{\dot{m}}_{\mathrm{GSHP}}+{\dot{m}}_{\mathrm{ASHP}})C_{\mathrm{w}}{T}_3+Q_3+Q_{\mathrm{case},3},\hfill \end{array}$$

(23)

$$\begin{array}{cc}\hfill M_4{C}_{\mathrm{w}}{\dot{T}}_4=& S_{\mathrm{CHP}}{Q}_{\mathrm{CHP}}/4+({\dot{m}}_{\mathrm{DHW}}+{\dot{m}}_{\mathrm{SH}})C_{\mathrm{w}}{T}_3+{\dot{m}}_{\mathrm{GSHP}}{C}_{\mathrm{w}}{T}_{\mathrm{GSHP},\mathrm{out}}+{\dot{m}}_{\mathrm{ASHP}}{C}_{\mathrm{w}}{T}_{\mathrm{ASHP},\mathrm{out}}\hfill \\ \hfill & -({\dot{m}}_{\mathrm{DHW}}+{\dot{m}}_{\mathrm{SH}}+{\dot{m}}_{\mathrm{GSHP}}+{\dot{m}}_{\mathrm{ASHP}})C_{\mathrm{w}}{T}_4+Q_4+Q_{\mathrm{case},4},\hfill \end{array}$$

(24)

$$M_5{C}_{\mathrm{w}}{\dot{T}}_5=S_{\mathrm{EH}}{Q}_{\mathrm{EH}}+{\dot{m}}_{\mathrm{DHW}}{C}_{\mathrm{w}}(T_4-T_5)+Q_5+Q_{\mathrm{case},5},$$

(25)

where the heat transfer through the vertical conduction of node i, represented as $Q_i$ (W) can be modeled as [89,116]:

$$Q_1=A_{\mathrm{HWST}}{k}_{\mathrm{w}}(T_2-T_1)/H_{\mathrm{HWST}},$$

(26)

$$Q_i=A_{\mathrm{HWST}}{k}_{\mathrm{w}}(T_{i+1}-T_i)/H_{\mathrm{HWST}}+A_{\mathrm{HWST}}{k}_{\mathrm{w}}(T_{i-1}-T_i)/H_{\mathrm{HWST}},i=2,3,4,$$

(27)

$$Q_5=A_{\mathrm{HWST}}{k}_{\mathrm{w}}(T_4-T_5)/H_{\mathrm{HWST}},$$

(28)

where $A_{\mathrm{HWST}}$ (m${}^2$) is the cross sectional area between the neighbouring nodes, $k_{\mathrm{w}}$ (W/(m K)) is the thermal conductivity of water and $H_{\mathrm{HWST}}$ (m) is the height of each layer. The heat transfer between water in node i and the inner surface of the tank $Q_{\mathrm{case},i}$ (W) can be calculated with Newton’s law of cooling:

$$Q_{\mathrm{case},i}=h_{\mathrm{case}}\left(\pi d_{\mathrm{HWST},\mathrm{in}}{H}_{\mathrm{HWST}}\right)\left(T_{\mathrm{case}}-T_i\right),$$

(29)

where the evolution of the case temperature $T_{\mathrm{case}}$ (${}^{\circ }\mathrm{C}$) is:

$$M_{\mathrm{case}}{C}_{\mathrm{case}}{\dot{T}}_{\mathrm{case}}=-\sum _{i=1}^5{Q}_{\mathrm{case},i}-Q_{\mathrm{HWST},\mathrm{L}},$$

(30)

where $Q_{\mathrm{HWST},\mathrm{L}}$ (W) is the heat loss of the tank calculated by:

$$Q_{\mathrm{HWST},\mathrm{L}}=U_{\mathrm{HWST},\mathrm{L}}\left(5\pi d_{\mathrm{HWST},\mathrm{out}}{H}_{\mathrm{HWST}}\right)(T_{\mathrm{case}}-T_{\mathrm{zone}}).$$

(31)

It is worth mentioning that considering the connection between HWST and other heating or cooling technologies is what allows us to build a multi-device framework combining several technologies.

3.7. Phase Change Materials

PCMs are becoming more popular and approachable for thermal storage due to their high storage density during the phase change process. The most common type of PCM is the solid–liquid one [118], which can be installed in an HWST [66,89], water-based heating/cooling system [119,120], ventilation system [121] or building structure [60,122]. Different kinds of PCM are chosen according to their physical characteristics, such as their heat convective coefficient and phase change temperature, for different application requirements. Water or air usually takes the role of the heat transfer fluid in the charging and discharging process.

The effective heat capacity method [60] and the enthalpy method [89] are two common ways to model the thermal storage process of PCM [123]. The charging performance of a PCM tank based on thermodynamic laws was studied in [66], whereas a two-dimensional PCM model in HWST based on the enthalpy method was studied in [89]. The model of the PCM-to-air heat exchanger can be found in [122]. An interesting comparison between the effective heat capacity method and enthalpy method can be found in [123]. It is shown that the two methods are both effective in predicting the charging or discharging time and air temperature, whereas the effective heat capacity method is more accurate in modeling the component temperatures of PCM. Thus, both of them can be used when studying the system level performance, but the effective heat capacity method is preferred when focusing on PCM behaviour.

In this study, it is assumed that the PCM has a unified temperature, demonstrating a high performance with no complex calculations [121]. The effective heat capacity method is expressed as [64]:

$$C_{\mathrm{PCM}}=\left\{\begin{array}{cc}{C}_{\mathrm{PCM},\mathrm{s}}& T_{\mathrm{PCM}}<T_{\mathrm{m}1}\\ C_{\mathrm{PCM},\mathrm{s}}+\frac{h_{\mathrm{PCM}}}{T_{\mathrm{m}2}-T_{\mathrm{m}1}}& T_{\mathrm{m}1}<T_{\mathrm{PCM}}<T_{\mathrm{m}2}\\ C_{\mathrm{PCM},\mathrm{l}}& T_{\mathrm{PCM}}>T_{\mathrm{m}2}\end{array}\right.,$$

(32)

$$M_{\mathrm{PCM}}{C}_{\mathrm{PCM}}{\dot{T}}_{\mathrm{PCM}}=Q_{\mathrm{PCM}},$$

(33)

where $C_{\mathrm{PCM},\mathrm{s}}$ (J/(kg K)) and $C_{\mathrm{PCM},\mathrm{l}}$ (J/(kg K)) are the specific heat of the PCM under a solid and liquid state, respectively, $h_{\mathrm{PCM}}$ (J/kg) is the latent heat of fusion, $T_{\mathrm{PCM}}$ (${}^{\circ }\mathrm{C}$) is the PCM temperature, $T_{\mathrm{m}1}$ (${}^{\circ }\mathrm{C}$) and $T_{\mathrm{m}2}$ (${}^{\circ }\mathrm{C}$) are the melting ranges and $Q_{\mathrm{PCM}}$ (W) is the heat transfer to PCM. The alternative enthalpy method is expressed as [89]:

$$h_{\mathrm{PCM}}=\left\{\begin{array}{cc}{C}_{\mathrm{PCM},\mathrm{s}}{T}_{\mathrm{PCM}}& T_{\mathrm{PCM}}<T_{\mathrm{m}1}\\ C_{\mathrm{PCM},\mathrm{s}}{T}_{\mathrm{m}1}+\frac{h_{\mathrm{sl}}\left(T_{\mathrm{PCM}}-T_{\mathrm{m}1}\right)}{T_{\mathrm{m}2}-T_{\mathrm{m}1}}& T_{\mathrm{m}1}\le T_{\mathrm{PCM}}\le T_{\mathrm{m}2}\\ C_{\mathrm{PCM},\mathrm{l}}\left(T_{\mathrm{PCM}}-T_{\mathrm{m}2}\right)+h_{\mathrm{sl}}+C_{\mathrm{PCM},\mathrm{s}}{T}_{\mathrm{m}1}& T_{\mathrm{PCM}}>T_{\mathrm{m}2}\end{array}\right.,$$

(34)

$${\rho }_{\mathrm{PCM}}{V}_{\mathrm{PCM}}{\dot{h}}_{\mathrm{PCM}}=Q_{\mathrm{PCM}},$$

(35)

where ${\rho }_{\mathrm{PCM}}$ (kg/m${}^3$), $V_{\mathrm{PCM}}$ (m${}^3$) and $h_{\mathrm{PCM}}$ (J/kg) are the density, volume and specific enthalpy of the PCM material. These two recommended models can support different application scenarios: it is necessary to couple the target technology (e.g., HWST) with the appropriate capacity (or enthalpy) coefficient resulting from the models above. This also allows us to combine different technologies into one framework.

4. Heating and Cooling Systems

According to the literature review, environmentally friendly heating and cooling energy technologies are inefficiently instructed in the state-of-the-art; therefore, they are taken into account in this section. They can be integrated to form the water-based energy system in buildings based on the water flow rate and/or temperature with the intention of the modeling process in a programming environment together with the control laws. This section concludes a review of the mathematical models of the environmentally friendly heating technologies (ASHP, GSHP), environmentally friendly cooling technologies (DEC, IEC, CC, AC, CCHP) and energy distribution systems (hydronic (HR) heating system, chilled beam, Radiant Cooling System (RCS)).

4.1. Heating Technologies

4.1.1. Air Source Heat Pump

ASHP transfers heat from the heat source to the heat sink depending on its mode, and this provides more potential for energy saving. Two types of ASHP are common nowadays [124]: one is the air-to-air heat pump, which transfers the heat extracted from outside air to inside air, and the other is the air-to-water heat pump, which transfers the heat extracted from ambient air to the HR heating system inside the building. The latter type is more convenient to integrate with the HWST and the water-based heating and cooling systems if they are available. Thus, this type of ASHP is considered for this study.

Heat pumps are usually modeled in terms of COP. The COP is defined as the ratio between the produced thermal power and consumed electric power. It is a widely used description for the efficiency of heat pumps and chillers. The COP of ASHP has been modeled in various ways, such as being a constant independent on the ambient temperature, supply water temperature and other factors [71], or a linear model of the ambient temperature $T_{\mathrm{amb}}$ [70]. The most accurate but complex model is a quadratic form, as studied in [125], taking into account the ambient temperature, normalized compressor power and returned water temperature to the condenser. A linear approximation of ${\mathrm{COP}}_{\mathrm{ASHP}}$ is used in this study [126]:

$${\mathrm{COP}}_{\mathrm{ASHP}}=6.1189+0.0676T_{\mathrm{amb}}-0.0632T_{\mathrm{ASHP},\mathrm{in}},$$

(36)

where $T_{\mathrm{ASHP},\mathrm{in}}$ (${}^{\circ }\mathrm{C}$) is the inlet water temperature. Then, the heating power $Q_{\mathrm{ASHP}}$ (W) and outlet temperature $T_{\mathrm{ASHP},\mathrm{out}}$ (${}^{\circ }\mathrm{C}$) can be calculated by:

$$Q_{\mathrm{ASHP}}={\mathrm{COP}}_{\mathrm{ASHP}}{W}_{\mathrm{ASHP}}={\dot{m}}_{\mathrm{ASHP}}{C}_{\mathrm{w}}(T_{\mathrm{ASHP},\mathrm{out}}-T_{\mathrm{ASHP},\mathrm{in}}),$$

(37)

where $W_{\mathrm{ASHP}}$ (W) is the power consumed by the compressor and ${\dot{m}}_{\mathrm{ASHP}}$ (kg/s) is the mass flow rate of water, which is usually assumed to be 0.5 kg/s. This model is a good trade-off between accuracy and complexity. It can be connected with other parts of the water-based heating systems and contributes to the whole framework.

4.1.2. Ground Source Heat Pump

GSHPs use ground temperature to warm or cool the circulating water in the system. Due to the fact that the ground temperature is more stable than the ambient air temperature, GSHP can usually achieve a better performance (i.e., higher COP) than ASHP.

The factors that influence the COP of a GSHP system include the type of ground heat exchangers (horizontal or vertical), soil condition, local climate, building load characteristics and system design and control optimization [65]. Again, the COP has been modeled as a constant [11], or dependent on meteorological information, such as solar radiation and the ambient temperature [127]. A field study conducted in hot summer and cold winter areas in China [128] indicated that the COP ranges from 1.95 to 4.35 depending on various outside conditions. A quadratic model for the COP of a solar-assisted GSHP was considered in [129], whereas the adverse influence of soil temperature evolution was quantified in [130], since the operation of GSHP caused variation in the soil temperature when the load was high. It is also shown that increasing the vertical length and the space between the exchangers can decrease the soil temperature variation. It is assumed that the soil temperature is constant in one month as a conclusion of the various proposed models, and that the quadratic model of ${\mathrm{COP}}_{\mathrm{GSHP}}$ based on the temperature difference between the soil and the outlet water $\Delta T_{\mathrm{GSHP}}$ is as follows [131]:

$${\mathrm{COP}}_{\mathrm{GSHP}}=8.77-0.150\Delta T_{\mathrm{GSHP}}+0.000734\Delta T_{\mathrm{GSHP}}^2,\mathrm{for}20\le \Delta T_{\mathrm{GSHP}}\le 60.$$

(38)

The operation of GSHP can be modeled similar to ASHP:

$$Q_{\mathrm{GSHP}}={\mathrm{COP}}_{\mathrm{GSHP}}{W}_{\mathrm{GSHP}}={\dot{m}}_{\mathrm{GSHP}}{C}_{\mathrm{w}}(T_{\mathrm{GSHP},\mathrm{out}}-T_{\mathrm{GSHP},\mathrm{in}}),$$

(39)

where $Q_{\mathrm{GSHP}}$ (W) and $W_{\mathrm{GSHP}}$ (W) are the heating and consumed compressor power, and ${\dot{m}}_{\mathrm{GSHP}}$ (kg/s) is the water mass flow rate, which is assumed to be 0.5 kg/s [12]. The modeling methods based on COP avoid the long calculations in the heat transfer models of ground heat exchanger; for example, in [132].

4.2. Cooling Technologies

Although novel cooling systems are continuously proposed as a combination of various technologies, exploring all of these combinations is outside the scope of this work. Therefore, a few popular technologies that are suitable for the aimed framework are carefully investigated.

4.2.1. Evaporative Cooling

Evaporative cooling is an efficient and sustainable way to cool the space in a hot and dry climate. It utilizes the evaporation process while the non-saturated water absorbs a large amount of heat from the air, increasing the air humidity and decreasing the air temperature simultaneously. Evaporative cooling is categorized as direct and indirect evaporative cooling (DEC and IEC) [133]. In DEC, water and air make contact directly, allowing the air to cool to wet bulb temperature, but the increased humidity of the outlet air may bring discomfort. In IEC, there is no direct contact between the water and air, allowing the air to cool to dew point temperature.

A mathematical model based on energy conservation, heat transfer and mass transfer for DEC was investigated in [134], calculating the total heat transfer with temperature difference and absolute humidity difference. This model is further simplified in [73] to improve the computational effectiveness:

$${\epsilon }_{\mathrm{DEC}}=1-\exp \left\{-\frac{h_{\mathrm{DEC}}{A}_{\mathrm{DEC}}}{{\dot{m}}_{\mathrm{DEC}}{C}_{\mathrm{hum}}}\right\},$$

(40)

where ${\epsilon }_{\mathrm{DEC}}$ (%) is the effectiveness of DEC, $h_{\mathrm{DEC}}$ (W/(m${}^2$ K)) is the convective heat transfer coefficient, $A_{\mathrm{DEC}}$ (m${}^2$) is the area of the heat transfer surface, ${\dot{m}}_{\mathrm{DEC}}$ (kg/s) is the mass flow rate of the air and $C_{\mathrm{hum}}$ (J/(kg K)) is the specific heat of the humid air, expressed as:

$$C_{\mathrm{hum}}=C_{\mathrm{air}}+H_{\mathrm{w}}{C}_{\mathrm{vap}},$$

(41)

where $C_{\mathrm{vap}}$ (J/(kg K)) is the specific heat of the vapour, and $H_{\mathrm{w}}$ (kg water/kg dry air) is the air humidity ratio of the inlet air (e.g., the ambient air). Typically, it is assumed that $h_{\mathrm{DEC}}$ = 35 W/(m${}^2$ K), ${\dot{m}}_{\mathrm{DEC}}$ = 1 kg/s (the heat transfer coefficient is based on the air velocity), $C_{\mathrm{air}}$=1007 J/(kg K) and $C_{\mathrm{vap}}$ = 1868 J/(kg K). Then, the outlet air temperature and total cooling power $Q_{\mathrm{DEC}}$ (W) are:

$$T_{\mathrm{DEC},\mathrm{out}}=T_{\mathrm{DEC},\mathrm{in}}-{\epsilon }_{\mathrm{DEC}}(T_{\mathrm{DEC},\mathrm{in}}-T_{\mathrm{wbt}}),$$

(42)

$$Q_{\mathrm{DEC}}={\dot{m}}_{\mathrm{DEC}}{C}_{\mathrm{air}}(T_{\mathrm{DEC},\mathrm{in}}-T_{\mathrm{DEC},\mathrm{out}}),$$

(43)

where $T_{\mathrm{DEC},\mathrm{in}}$ (${}^{\circ }\mathrm{C}$) and $T_{\mathrm{DEC},\mathrm{out}}$ (${}^{\circ }\mathrm{C}$) are the temperatures of the inlet and outlet air respectively, and $T_{\mathrm{wbt}}$ (${}^{\circ }\mathrm{C}$) is the wet bulb temperature of the air.

The Number of Transfer Unit (NTU) method, also referred to as the $\epsilon$-NTU method, is used to deal with the situation where the outlet temperature is unknown, such as evaporative cooling systems and chilled beams. An IEC model based on a modification of this method [74] was proven to significantly reduce the computational complexity compared to the differential equation models based on physical calculations. This method assumed a linear relationship between the temperature of saturated air in the water film $T_{\mathrm{sat}}$ (${}^{\circ }\mathrm{C}$) and the enthalpy h (J/kg) within a small temperature range:

$$h\left(T_{\mathrm{sat}}\right)=a{T}_{\mathrm{sat}}+b,$$

(44)

where the values of a and b can be fitted with the data in [135]. Then, the cooling power $Q_{\mathrm{IEC}}$ (W) is expressed as:

$$Q_{\mathrm{IEC}}={\epsilon }_{\mathrm{IEC}}^{*}{C}_{\mathrm{IEC},\min }\left(h\left(T_{\mathrm{dry},\mathrm{in}}\right)-h_{\mathrm{wet},\mathrm{in}}\right),$$

(45)

where $T_{\mathrm{dry},\mathrm{in}}$ (${}^{\circ }\mathrm{C}$) is the inlet air temperature to the dry passage, $h_{\mathrm{wet},\mathrm{in}}$ (J/kg) is the enthalpy of inlet air to the wet passage, ${\epsilon }_{\mathrm{IEC}}^{*}$ (%) is the modified effectiveness and $C_{\mathrm{IEC},\min }$ is calculated based on the modified cold and hot fluid heat capacity:

$$C_{\mathrm{c}}^{*}={\dot{m}}_{\mathrm{wet}},$$

(46)

$$C_{\mathrm{h}}^{*}={\dot{m}}_{\mathrm{dry}}{C}_{\mathrm{air}}/a,$$

(47)

$$C_{\mathrm{IEC},\min }=\min \left(C_{\mathrm{c}}^{*},C_{\mathrm{h}}^{*}\right),$$

(48)

where ${\dot{m}}_{\mathrm{wet}}$ (kg/s) and ${\dot{m}}_{\mathrm{dry}}$ (kg/s) are the air mass flow rate in the wet and dry passage. ${\epsilon }_{\mathrm{IEC}}^{*}$ is calculated based on a function of modified NTU values. For a parallel flow IEC, the function is:

$${\epsilon }_{\mathrm{IEC}}^{*}=\frac{1-\exp \left[-{\mathrm{NTU}}^{*}\left(1+C_{\mathrm{IEC},r}\right)\right]}{1+C_{\mathrm{IEC},r}}.$$

(49)

For a counter flow IEC, the function is:

$${\epsilon }_{\mathrm{IEC}}^{*}=\frac{1-\exp \left[-{\mathrm{NTU}}^{*}\left(1-C_{\mathrm{IEC},\mathrm{r}}\right)\right]}{1-C_{\mathrm{IEC},\mathrm{r}}\exp \left[-{\mathrm{NTU}}^{*}\left(1-C_{\mathrm{IEC},\mathrm{r}}\right)\right]}\mathrm{where}{C}_{\mathrm{IEC},\mathrm{r}}<1,$$

(50)

where the modified NTU values ${\mathrm{NTU}}^{*}$ and capacity ratio $C_{\mathrm{IEC},\mathrm{r}}$ are calculated by:

$${\mathrm{NTU}}^{*}=\frac{U^{*}{A}_{\mathrm{IEC}}}{C_{\mathrm{IEC},\min }}\mathrm{where}{U}^{*}=\frac{1}{a(\frac{1}{h_{\mathrm{IEC}}}+\frac{{\delta }_{\mathrm{IEC}}}{k_{\mathrm{IEC}}})+\frac{1}{{\beta }_{\mathrm{IEC}}}},$$

(51)

$$C_{\mathrm{IEC},\max }=\max \left(C_{\mathrm{c}}^{*},C_{\mathrm{h}}^{*}\right),$$

(52)

$$C_{\mathrm{IEC},\mathrm{r}}=C_{\mathrm{IEC},\min }/C_{\mathrm{IEC},\max },$$

(53)

where $h_{\mathrm{IEC}}$ (W/(m${}^2$ K)) is the convective heat transfer coefficient in the dry passage, ${\delta }_{\mathrm{IEC}}$ (m) and $k_{\mathrm{IEC}}$ (W/(m K)) are the combined thickness and thermal conductivity of the water film and thin wall between the dry and wet passages and ${\beta }_{\mathrm{IEC}}$ (kg water/kg dry air) is the mass transfer coefficient in the wet passage. The typical parameter settings for a 0.25 m${}^2$ IEC are ${\dot{m}}_{\mathrm{wet}}$ = 0.00098 kg/s, ${\dot{m}}_{\mathrm{dry}}$ = 0.0014 kg/s and ${\delta }_{\mathrm{IEC}}=0.0015$ m. These two models are recommended for DEC and IEC since they can be used for different configurations without considering the dynamic energy transfer process.

4.2.2. Chiller

The cooling process of a chiller is similar to the inverse of the heating process of a heat pump. There are two kinds of chillers: absorption and compression chillers (ACs and CCs) [81]. The first type uses thermal heat to circulate the refrigerant around the system, so that it is profitable in the situation where heat can be available at low prices, e.g., the waste heat from the boiler. The second type uses a compressor to circulate the refrigerant. It is suitable for the situation where electricity is cheap and heat is expensive to achieve. COP is the common performance factor for chillers. Similar to the COP of heat pumps, the COP of chillers also can be modeled as a constant [136], a linear function of system states [137,138], a function of the part load ratio [139], a bi-quadratic function [80], polynomials [80] or an exponential function [79].

Four modeling methods for ACs, including a white-box, a grey-box and two black-box models, are compared in [79], indicating that the empirical models can offer a simple way to calculate the cooling capacity and COP without knowing the internal working fluid temperature. Six black-box models based on water-side data for CCs were compared in [80], showing that the bi-quadratic model and polynomial model had a better performance. An accurate linear COP approximation method of chillers was proposed in [137]: the inputs include the inlet and outlet water temperatures, source and sink temperatures and the cooling load. To avoid a large number of inputs that need extra measurements, COP models based on part load ratios are adopted in this work. The model for AC is presented in [41] as follows:

$${\mathrm{COP}}_{\mathrm{AC}}=-0.6181{\eta }_{\mathrm{AC}}^2+0.8669{\eta }_{\mathrm{AC}}+0.4724.$$

(54)

The model for CC is presented in [83] as follows:

$${\mathrm{COP}}_{\mathrm{CC}}={\epsilon }_{\mathrm{CC}}(2.886·{10}^{-3}{P}_{\mathrm{CC}}^2+0.0293P_{\mathrm{CC}}+4.711),$$

(55)

$${\epsilon }_{\mathrm{CC}}=-0.569{\eta }_{\mathrm{CC}}^3-0.258{\eta }_{\mathrm{CC}}^2+1.520{\eta }_{\mathrm{CC}}+0.321,$$

(56)

where ${\eta }_{\mathrm{AC}}$ (%) and ${\eta }_{\mathrm{CC}}$ (%) are the part load ratios of AC and CC, and these two values usually range from 0.2 to 1 due to the minimum engine speed. ${\epsilon }_{\mathrm{CC}}$ (%) is the relative efficiency of correcting the COP of the chiller under different part load ratios, and $P_{\mathrm{CC}}$ (W) is the capacity of CC. Then, the cooling power and returned chilled water temperature can be calculated by:

$$Q_{\mathrm{AC}}={\mathrm{COP}}_{\mathrm{AC}}{W}_{\mathrm{AC}}={\dot{m}}_{\mathrm{AC}}{C}_{\mathrm{w}}(T_{\mathrm{AC},\mathrm{out}}-T_{\mathrm{AC},\mathrm{in}}),$$

(57)

$$Q_{\mathrm{CC}}={\mathrm{COP}}_{\mathrm{CC}}{W}_{\mathrm{CC}}={\dot{m}}_{\mathrm{CC}}{C}_{\mathrm{w}}(T_{\mathrm{CC},\mathrm{out}}-T_{\mathrm{CC},\mathrm{in}}).$$

(58)

The CCHP system is usually a combination of CHP and chillers. The generated energy by CHP in terms of electricity or heat can be supplied to chillers to provide cooling energy [44,140]. The connections and energy flows in a CCHP system are shown in Figure 3. The power generation unit produces electricity with fuel inputs (e.g., gas), and the recovered heat is collected. The generated electric power $P_{\mathrm{CHP}}$ (W) and thermal power $Q_{\mathrm{CHP}}$ (W) can be calculated with the CHP model presented in Section 3.4. Then, the generated electricity and heat, possibly along with other supplements from energy storage devices or the grid, support the CCs or ACs. Thus, the final generated cooling power of CCHP can be calculated with the chiller model presented in Section 4.2.2 with the electricity or thermal inputs. The performance of CCHP is evaluated by the efficiency of the power generation unit and the COP of the chillers [41].

4.3. Energy Distribution Systems

The supplied energy is distributed into each zone by the hydronic distribution systems used in this research. The models of auxiliary systems (e.g, water pumps) placed between energy systems and distribution systems are not studied in this work.

4.3.1. Hydronic Heating Systems

HR heating is a popular heat distributed system as it is compatible with several energy sources, such as solar energy and heat pumps. Combining conduction and radiant technology, heat emitters, such as wall-mounted radiators and underfloor coils, emit radiant heat to provide heat [141]. The output power $Q_{\mathrm{HR}}$ (W) of a HR radiator can be modeled as follows [142]:

$$Q_{\mathrm{HR}}=K_{\mathrm{HR}}(T_{\mathrm{HR},\mathrm{w}}-T_{\mathrm{zone}}),$$

(59)

$$M_{\mathrm{HR},\mathrm{w}}{C}_{\mathrm{w}}{\dot{T}}_{\mathrm{HR},\mathrm{w}}={\dot{m}}_{\mathrm{HR}}{A}_{\mathrm{HR}}{C}_{\mathrm{w}}{T}_{\mathrm{HR},\mathrm{in}}-{\dot{m}}_{\mathrm{HR}}{A}_{\mathrm{HR}}{C}_{\mathrm{w}}{T}_{\mathrm{HR},\mathrm{out}}-Q_{\mathrm{HR}},$$

(60)

where $M_{\mathrm{HR},\mathrm{w}}$ (kg) and $T_{\mathrm{HR},\mathrm{w}}$ (${}^{\circ }\mathrm{C}$) are the mass and averaged temperature of water inside the HR radiator, $A_{\mathrm{HR}}$ (m${}^2$) is the area of the radiator and $K_{\mathrm{HR}}$ (W/K) is the heat transfer coefficient that depends on the temperature difference. For simplicity, it is assumed that the radiator surface and inside water have the same temperature, so that $K_{\mathrm{HR}}$ can be calculated as:

$$K_{\mathrm{HR}}=h_{\mathrm{HR}}{A}_{\mathrm{HR}},$$

(61)

where $h_{\mathrm{HR}}$ (W/(m${}^2$ K)) is the convection coefficient, and $A_{\mathrm{HR}}$ (m${}^2$) is the area of the radiator. The typical parameter settings of a low temperature HR radiator are $h_{\mathrm{HR}}=4.84$ W/(m${}^2$ K) and ${\dot{m}}_{\mathrm{HR}}=0.04$ kg/(m${}^2$ s) [143]. This model is based on an energy balance equation and can be connected with the models of thermal storage or heating technologies presented in this work.

4.3.2. Chilled Beam

The chilled beam is a cooling distribution technology utilized in a wide range of types of buildings, such as industries, commercial buildings, households and even transportation vehicles. There are mainly two kinds of chilled beams: Passive Chilled Beams (PCBs) and Active Chilled Beams (ACBs) [144]. The schematic representations of them are shown in Figure 4 and Figure 5.

The PCB provides cooling power for buildings through natural convection. The rising warm air loses heat to the coils with chilled water and falls to cool the room. The cooling capacity of PCB is dependent on the temperature difference (typically 8–9 ${}^{\circ }\mathrm{C}$) between the chilled water and indoor air, beam configuration and height of the roof [145]. Typically, the chilled water at 14 ${}^{\circ }\mathrm{C}$ to 18 ${}^{\circ }\mathrm{C}$ flowing through the coils at speed ranges from 0.02 kg/s to 0.2 kg/s and provides 20–40 W/(m K) cooling power [146].

Analytical PCB models based on CFD analysis can accurately predict the indoor airflow, air temperature and contaminant concentration [147], but they are time consuming. Alternatively, the simplified mathematical models that are easily implemented can be developed by a regression method. For example, a multi-variable regression model with the data from the CFD analysis for the coupled displacement-ventilation and PCB system was investigated in [148]. Another work based on similar methods with the data from experimental studies predicted the cooling capacity $Q_{\mathrm{PCB}}$ (W) based on part load ratios [90], expressing the models for a single chilled beam as:

$$\frac{Q_{\mathrm{PCB}}}{Q_{\mathrm{PCB},\mathrm{rated}}}=0.712\times {\left(\frac{{\dot{m}}_{\mathrm{PCB}}}{{\dot{m}}_{\mathrm{PCB},\mathrm{rated}}}\right)}^{0.707}\times {\left(0.142+\frac{T_{\mathrm{zone}}-T_{\mathrm{PCB},\mathrm{in}}}{T_{\mathrm{rated}}-T_{\mathrm{PCB},\mathrm{in},\mathrm{rated}}}\right)}^{1.714},$$

(62)

and the model for multiple chilled beams in a large open office is:

$$\frac{Q_{\mathrm{PCB}}}{Q_{\mathrm{PCB},\mathrm{rated}}}=0.612\times {\left(\frac{{\dot{m}}_{\mathrm{PCB}}}{{\dot{m}}_{\mathrm{PCB},\mathrm{rated}}}\right)}^{0.464}\times {\left(0.132+\frac{T_{\mathrm{zone}}-T_{\mathrm{PCB},\mathrm{in}}}{T_{\mathrm{rated}}-T_{\mathrm{PCB},\mathrm{in},\mathrm{rated}}}\right)}^{0.866},$$

(63)

where the subscript $\mathrm{rated}$ means the rating case, and the typical values are: $Q_{\mathrm{PCB},\mathrm{rated}}=703$ W, ${\dot{m}}_{\mathrm{PCB},\mathrm{rated}}=0.08$ kg/s and $T_{\mathrm{rated}}=22.5{}^{\circ }\mathrm{C}$. These formulations simplified the calculation by avoiding the physical details in heat transfer and distribution, and have been used in a building level analysis in [149].

The ACB provides cooling power for the building through forced convection. The indoor air is inducted to flow across the coil with chilled water. Typically, the induction ratio is 1/4 and the cooling power is 100–200 W/(m K) [146]. The NTU model for ACB presented in [150] is proven to be accurate and easy to implement for modeling the cooling capacity and supplied air temperature compared with other accurate analytical models; for example, in [151], it expresses the total cooling power $Q_{\mathrm{ACB}}$ (W) by:

$$Q_{\mathrm{ACB}}=Q_{\mathrm{pri}}+Q_{\mathrm{w}},$$

(64)

where $Q_{\mathrm{pri}}$ (W) is the cooling power from the primary air:

$$Q_{\mathrm{pri}}={\dot{m}}_{\mathrm{pri}}{C}_{\mathrm{air}}(T_{\mathrm{pri}}-T_{\mathrm{exh}}),$$

(65)

where ${\dot{m}}_{\mathrm{pri}}$ (kg/s) is the mass flow rate of the primary air flow, which is typically equal to 0.03 kg/s, and $T_{\mathrm{pri}}$ (${}^{\circ }\mathrm{C}$) and $T_{\mathrm{exh}}$ (${}^{\circ }\mathrm{C}$) are the temperatures of the primary air and exhausted air. $Q_{\mathrm{w}}$ (W) is the cooling power from chilled water, and is calculated by the NTU method:

$$Q_{\mathrm{w}}={\epsilon }_{\mathrm{ACB}}{C}_{\mathrm{ACB},\min }(T_{\mathrm{zone}}-T_{\mathrm{ACB},\mathrm{in}}),$$

(66)

where $T_{\mathrm{ACB},\mathrm{in}}$ (${}^{\circ }\mathrm{C}$) is the inlet temperature of the chilled water, and $C_{\mathrm{ACB},\min }$ (J/(kg K)) is the minimum heat capacity according to:

$$C_{\mathrm{ACB},\min }=\min ({\dot{m}}_{\sec }{C}_{\mathrm{air}},{\dot{m}}_{\mathrm{ACB}}{C}_{\mathrm{w}}),$$

(67)

where ${\dot{m}}_{\sec }$ (kg/s) refers to the mass flow rate of secondary air (i.e., induced air), which is the product of the induction ratio, and ${\dot{m}}_{\mathrm{pri}}$, ${\dot{m}}_{\mathrm{ACB}}$ (kg/s) is the mass flow rate of the chilled water, which is typically equal to 0.04 kg/s. ${\epsilon }_{\mathrm{ACB}}$ (%) is the effectiveness parameter, which is a function of NTU, for the common application cases (NTU is less than 3) of cross-flow ACB, and an accurate approximation of ${\epsilon }_{\mathrm{ACB}}$ is:

$${\epsilon }_{\mathrm{ACB}}=\left(\frac{\mathrm{NTU}}{1+1.1238·\mathrm{NTU}}+e^{-\mathrm{NTU}}-1\right)·C_{\mathrm{ACB},r}+1-e^{-\mathrm{NTU}},$$

(68)

where $C_{\mathrm{ACB},r}$ is calculated with:

$$C_{\mathrm{ACB},\max }=\max (m_{\sec }{C}_{\mathrm{air}},m_{\mathrm{ACB}}{C}_{\mathrm{w}}),$$

(69)

$$C_{\mathrm{ACB},r}=C_{\mathrm{ACB},\min }/C_{\mathrm{ACB},\max },$$

(70)

and NTU is defined as:

$$\mathrm{NTU}=\frac{h_{\mathrm{ACB}}{A}_{\mathrm{ACB}}}{C_{\mathrm{ACB},\min }},$$

(71)

where $h_{\mathrm{ACB}}$ (W/(m${}^2$ K)) is the overall heat transfer coefficient of the coil and $A_{\mathrm{ACB}}$ (m${}^2$) refers to the heat transfer area. The NTU model is recommended for the aimed framework since it shows a good correspondence with measured data without extensive calibrations and long calculations.

4.3.3. Radiant Cooling System

Another popular terminal cooling technology is the RCS. It is a broad concept that refers to a system where radiant heat transfer covers over half of the total heat transfer. The embedded surface cooling system, thermally activated building system and radiant cooling panels all belong to this concept [152]. It has been increasingly applied in a wide range of types of buildings due to its lower energy consumption and better thermal comfort compared to traditional air conditioners [153]. A simplified model of RCS was proposed in [154]: assuming that all masses in the room have temperatures close to the zone air and that the slab temperature is uniform, then the cooling power $Q_{\mathrm{RCS}}$ (W) is calculated as follows:

$$Q_{\mathrm{RCS}}=U_{\mathrm{RCS}}{A}_{\mathrm{RCS}}(T_{\mathrm{zone}}-T_{\mathrm{RCS}}),$$

(72)

where $U_{\mathrm{RCS}}$ (W/(m${}^2$ K)) is the heat transfer coefficient expressed by:

$$U_{\mathrm{RCS}}=k{|T_{\mathrm{zone}}-T_{\mathrm{RCS}}|}^{0.09},$$

(73)

where $T_{\mathrm{RCS}}$ (${}^{\circ }\mathrm{C}$) is the radiant surface temperature, and k is a constant depends on the emissivities of the panel. The values can be found in [155]; one typical value is 8. The dynamics of the slab temperature $T_{\mathrm{RCS}}$ is:

$$M_{\mathrm{RCS}}{C}_{\mathrm{RCS}}{\dot{T}}_{\mathrm{RCS}}={\dot{m}}_{\mathrm{RCS}}{C}_{\mathrm{w}}{T}_{\mathrm{RCS},\mathrm{in}}-{\dot{m}}_{\mathrm{RCS}}{C}_{\mathrm{w}}{T}_{\mathrm{RCS},\mathrm{out}}+Q_{\mathrm{RCS}}+{\alpha }_{\mathrm{int}}{Q}_{\mathrm{int}},$$

(74)

where $M_{\mathrm{RCS}}$ (kg) and $C_{\mathrm{RCS}}$ (J/(kg K)) are the mass and specific heat of the radiant panel, ${\dot{m}}_{\mathrm{RCS}}$ (kg/s) is the water mass flow rate, $T_{\mathrm{RCS},\mathrm{in}}$ (${}^{\circ }\mathrm{C}$) and $T_{\mathrm{RCS},\mathrm{out}}$ (${}^{\circ }\mathrm{C}$) are the inlet and outlet temperatures of the panel and ${\alpha }_{\mathrm{int}}$ (%) is the fraction of the internal heat gain absorbed by the radiant panel. A typical value for ${\dot{m}}_{\mathrm{RCS}}$ is 0.08 kg/s for every 100 m pipe. The model is validated under weekly prediction and the complexity allows for its integration with a real-time MPC algorithm [154].

5. Energy Management System

In the previous sections, building physics and energy systems were discussed to create a framework for buildings equipped with several energy systems. In this section, the picture of the whole framework is completed by discussing how to control these technologies, i.e., EMS is investigated. Some popular control strategies are reviewed in the following stages of this work to show how to operate the technologies based on these models. This investigation reviews EMS approaches to categorize them in terms of their visions and usability.

A solution to enable EMS is to design a combination of computer-aided strategies to monitor and control building energy loads [156]. This system aims to improve the performance efficiency of the installed technologies, including the HVAC devices, onsite generation and storage systems, by operating them correctly. In this way, EMS allows buildings to achieve both their own objectives and those of grid suppliers based on the prediction of the energy generation and demand [157]. Two strategies are most frequently used to design the EMS: Rule-Based Control (RBC) and Model Predictive Control (MPC), respectively. This section examines the formulation of control problems and the concepts, approaches and comparison of RBC and MPC.

5.1. Rule-Based Control

RBC is a series of predefined decisions based on “IF-THEN” commands and real-time monitoring. The “IF” conditions of these commands often include some threshold values dependent on the difference between a reference value and the measured values (e.g., limited price and hourly energy prices). RBC is adjusted by determining the operation of energy systems (e.g., on/off of HVAC devices) based on a rule set. The control objectives are set to reduce the energy cost and energy consumption, for example.

RBC is most frequently used for cost reduction based on real-time energy prices. A multi-level RBC algorithm for a solar-assisted GSHP heating system was investigated in [158], showing that the EMS could select different energy technologies under the fluctuating real-time electricity price. Three different RBC approaches as examples, called momentary, backwards-looking and predictive control algorithms, were proposed and compared in [56] to control a GSHP heating system. It was shown that predictive control outperforms the others. The control rules in the three approaches were based on the comparison between the current Hourly Energy or Electricity Price (HEP) and a constant limiting price, the comparison between the current HEP and the median of the past HEP and the trend of future HEP, respectively. The typical methods to generate the price trend include a blocking-maximum subarray method, sliding-maximum subarray method and moving average [21]. Different definitions of price thresholds for RBC based on the location and length of the time window were discussed in [159]. The price thresholds can be calculated with either the bottom and upper percentiles, or deviations from the average prices of the time window. In addition to controlling the heating or cooling systems as discussed above, RBC was also employed for controlling the operation of multi electricity appliances and the storage system [160], the setpoint temperature of swimming halls [161], the ${\mathrm{CO}}_2$ concentration of educational office buildings [162] and the integrated operation of HVAC and blinds in office buildings [163] with similar rules. Occupancy patterns can also be taken into consideration when setting the rules, with the set-back and pre-cooling modes defined based on the occupancy profile [164].

Although some of the above studies showed an over 10% cost reduction under real-time electricity prices, there are still plenty of situations that do not provide real-time energy prices, leading to a lower economic potential by load shifting. As a result, RBC is also implemented to decrease the energy consumption. A rule set to preheat or switch off the equipment based on the historical data was defined in [165], and another one to control the operations of smart windows to minimize the solar flux and guarantee the illumination in cooling seasons was studied in [166]. A fast chiller power control was employed in [167] with the rule set to estimate the chiller power reduction and to reset the air temperature setpoint and chilled water distribution. During transition seasons, RBC can be implemented by switching the HVAC system between three modes, called heating, passive cooling and active cooling, dependent on the thresholds of the ambient temperature and HWST temperature per hour [168], or activating the heating devices based on predefined schedules, determining the setpoint temperature based on a linear function of the ambient temperature and switching between two modes (winter and summer modes) based on the ambient temperature in the past 36 h [169]. The Rete algorithm was used to accelerate the rule processing when a large number of nodes are considered [170].

In addition, other objectives that RBC can deal with include maximizing the use of renewable energy by defining sequential priorities for a series of actions [171], reducing carbon emissions by determining the high-carbon threshold and low carbon threshold based on a sensitivity analysis [172], decreasing the peak-hour energy use by the pre-heating of the HWST and the building thermal mass [172] or charging and discharging the battery [173].

To conclude, RBC can either be as simple as a few rules or as complicated as rule sets, but all of them are in the form of straightforward rules or fundamental commands based on the monitoring of the system status and on predefined conditions. RBC is totally based on the current state, so less information, such as weather forecasting and load prediction, are needed. The extreme simplicity and easiness of implementation allow for its integration with analytical models for buildings and technologies [174]. Even if it depends on rules, it is useful in single buildings with few distributed generation technologies [160]. It is noted that once objective functions and scenarios are simple—for example, single zone, coherent objectives, or a limited number of design variables—RBC is able to generate results comparable to advanced strategies, such as MPC, because there are not sufficiently flexible for optimization [166].

However, the simplicity of RBC also brings some challenges. Firstly, there is no objective function in RBC formulation, so it is often used to deal with a single objective, and fails to deal with the trade-off between conflicting multi objectives (e.g., energy cost and thermal comfort), although, sometimes, some coherent objectives can be achieved simultaneously (e.g., reducing the energy cost and peak load [160]). Secondly, the different time constants of the hybrid production and emission systems and the randomness of the weather and occupancy bring challenges to tune the rules sets, which may degrade the performance of RBC, especially in large-scale non-residential buildings [168]. Thirdly, RBC cannot deal with the time delay due to the thermal inertia [175]: it tends to be out of operation and leads to energy waste when heating and cooling are needed simultaneously in transition periods. Last but not least, RBC is unavailable for providing optimal solutions, and there is still the significant saving potential of RBC results compared with that of advanced control strategies, especially in transient seasons [169]. Therefore, a more superior control is highly regarded, and MPC is one of the promising candidates. RBC can be used as a benchmark for MPC [173].

5.2. Model Predictive Control

An MPC formulation consists of a cost function to represent the control objective, several constraints for the system states and inputs, a predictive model (typically a linear model, since a nonlinear model may make MPC computationally cumbersome [28,34]) to mathematically predict the future system state, a finite predictive horizon and an optimization solver to generate control actions at each time step. The general form of MPC (in continuous time—discrete time implementation follows accordingly after discretization) is [34]:

$$\begin{array}{cc}\min \hfill & J={\int }_{\tau =t}^{t+N}f(x\left(\tau \right),u\left(\tau \right),y\left(\tau \right))d\tau \hfill \\ \mathrm{s}.\mathrm{t}.\hfill & x\in \mathcal{X}\hfill \\ \hfill & u\in \mathcal{U}\hfill \\ \hfill & \dot{x}=g(x,u,w)\hfill \\ \hfill & y=h\left(x\right)\hfill \end{array}$$

(75)

where $x$ is the system state, typically regarded as node temperatures, $u$ is the system input, $w$ is the system disturbance based on weather data, $\mathcal{X}$ and $\mathcal{U}$ are the constraints sets for the states and inputs, J is the cost function and N is the prediction horizon. The inputs of this framework are the current state $x$, the system model g and disturbances $w$. Every iteration generates an optimal input sequence of the input $u$ by solving the optimization problem, but only the first step result is implemented. Then, the prediction horizon recedes step by step and generates result for the new slot.

5.2.1. MPC Formulation

The control objectives of MPC are often expressed as penalties in the cost function J. The most frequently considered penalties in the objective functions include: (1) thermal comfort cost by penalizing the difference between the predictive and reference values of the Predicted Mean Vote (PMV) index [176], the deviation between the predictive and setpoint temperature [177] or the energy not supplied compared with the energy demand [178]; (2) the energy bill of the building [179]; (3) the energy consumption of the installed technologies [180]; (4) the carbon emissions calculated by the product of the carbon intensity and consumption of power [181]; (5) the nonrenewable energy proportion in the total load [7]. The MPC formulation is often deployed to find the optimal trade-off of the probable competitive objectives.

A variety of constraints for the state $\mathcal{X}$ and the input $\mathcal{U}$ can be handled by MPC. The constraints on the appliance operation power, thermal comfort and starting time of shiftable loads are typically represented as inequalities, whereas the constraints on building dynamics and the storage state are often represented as equalities. Slack variables, also called soft-constraints, can be included in the inequality constraints to avoid infeasibility [182,183].

The predictive model is regarded as the cornerstone of MPC [184]. An oversimplified predictive model may lead to extra tuning work and make the derived EMS unacceptable in real-world buildings [185], whereas the nonlinearity of an elaborate model may make the optimization step more complicated and the global optimum unavailable due to the non-convexity [186]. Thus, the predictive model should be designed carefully.

There are three kinds of modeling methods, which are called white-box (physics-based), black-box (data-driven) and grey-box (hybrid) methods, respectively. The white-box models describe the building energy performance based on the physical calculations, which need an exhaustive understanding of the technical characteristics [186]. For example, a high-fidelity white-box model is implemented in a large scale office building with a hybrid production system, and the optimization formulation is solved by decoupling the nonlinear problem to several smaller subproblems [168]. Sometimes, the white-box models refer to the well-established models in building simulation tools, such as TRNSYS and EnergyPlus [184].

The black-box models estimate the building energy performance dependent on numerical evaluations and statistical data. Two methods, called the autoregressive moving average with extra input model identification and subspace identification to generate linear black-box models, were compared in [187], demonstrating that both methods shows a high accuracy but that the latter is faster to implement due to fewer parameters. Alternatively, another work trained and refined the model by subspace system identification and a prediction error method based on excitation data during two weeks [188]. Nonlinear black-box models based on machine learning algorithms are also developed, including a dynamic artificial neural network to predict the building dynamics and thermal comfort for an office [189], regression trees and random forests to predict the room temperature and energy usage for a multi-zone commercial building [190] and genetic programming to predict the room temperature for a building with a conventional HR radiator supplied by a boiler [177].

The grey-box models are combinations of white-box and black-box models. They rely on a simplified relationship of a physical process, and the parameters are regressed based on measured data and have been applied to model single energy technology, such as condensing boilers [191], thermal zone and HVAC systems [192] and building physics [193]. A high order grey-box implementation can be found in [169]: a 7R-4C model was implemented for an office building with years of measurement data, and is able to be solved by the built-in solver in MATLAB.

5.2.2. Solution Techniques

The solution technique refers to the algorithm or strategy to solve the MPC formulation established in the previous section, which is in the form of an optimization problem. There are usually three kinds of solution techniques for these problems—implicit MPC, explicit MPC and approximate MPC—among which, the implicit way is the most frequently used one [194].

The idea of implicit MPC is to solve the optimization problem directly by the online searching of the control signal sequence over the predefined horizon. The computational complexity of this method depends on the MPC formulation discussed above (type of model, type of constraints, length of the horizon). Classical techniques, such as linear programming [7], quadratic programming [176], mixed integer programming [195] and dynamic programming [196], are the most notable methods for implicit MPC. Sometimes, hybrid optimization methods have been used to deal with nonlinearities, such as the combination of exhaustive search methods and subsequent quadratic programming used in [197].

The explicit MPC is developed to calculate the feasible values of the optimizer and to express the control variable $u$ explicitly as a function of the parameters. An example of its implementation is the scenario-based explicit MPC to decrease the energy use and to maintain the thermal and ${\mathrm{CO}}_2$ comfort level proposed in [198]: the control inputs are expressed as a function of the system states. However, the number of states is usually restricted to around 10 due to the computational complexity in calculating the explicit expression. This inherent weakness limits the use of explicit MPC only in simplified building models.

The approximate MPC trains a machine learning model to learn the MPC formulation and expresses $u$ as an implicit form of the parameters, predicting the control commands based on the same inputs of MPC without solving the optimization problem. One of the earliest implementations of this method in buildings is the one in [199], executing offline calculation integrated with building simulation tools. Two multivariate regression algorithms, called deep time-delay neural networks and regression trees, are studied in [200] to mimic the MPC controller, where the versatile framework is applicable without advanced software libraries. Another example is the recurrent neural network with a nonlinear autoregressive network with an exogenous input structure investigated in [189], demonstrating that approximate MPC can achieve most benefits of MPC while saving a large amount of computation time.

Although the building sectors can significantly benefit from the utilization of MPC, higher hardware and software levels are required [200,201], making it costly to implement. As a result of the complex calculation, most of the existing studies are carried out for linear predictive models, which restricts its performance and practical use [176]. MPC does not always dominate RBC. The comparison shows that MPC is more useful compared to simpler controllers, such as RBC, only when competitive objectives, multi-zone interactions and more design variables are considered [166].

Comprising the problem formulation and solution technique, MPC provides significant energy saving potential and performance enhancement, especially for complex building energy systems [199]. The objective function of MPC makes it possible to accommodate competitive control objectives; for example, to search for the trade-off between energy consumption and comfort [202]. Meanwhile, its ability to control multiple technologies, integrate the constraints and consider future predictions leads to a better performance compared to classical RBC strategies [178]. The predictive manner of MPC allows us to handle the disturbance with a robust or stochastic controller [31], and the receding horizon helps to avoid error accumulation [197]. These characteristics provide MPC with the ability to achieve 20% average cost savings and a 27% emission reduction [34].

To conclude, MPC is more attractive to handle complex objective functions, multiple zones and multiple technologies. The mathematical models presented in this work can be used to set up the optimization formulations according to the configurations of building and control objectives.

6. Discussion and Conclusions

Nowadays, new technologies, such as renewable energy resources and environmentally friendly technologies with a high efficiency, have been increasingly used to deal with environmental issues; for example, from the buildings’ side. EMSs provide effective ways to integrate and operate all of these technologies. They can improve the building performance significantly with well-designed laws or optimization methods. However, most of the complex EMSs studied in the literature still depend on the separate processes to model and control buildings in building simulation tools and programming environments, respectively, which brings challenges in setting up a real-time framework. No assemblable and computationally efficient mathematical models for popular energy systems have been investigated and applied at the building level. A coherent framework to integrate the applicable mathematical models of energy systems and control strategies is still an open problem.

This paper proposed the integration of mathematical models of building physics and energy systems in an integrative fashion. All of the technologies taken into account in this study proved to be beneficial in the literature review. The proposed models cover most of the popular technologies in buildings and avoid unnecessary details to make it simple when keeping the accuracy. They are categorized into energy generation, storage, consumption and distribution technologies to protect the buildings’ demand. The connections of these technologies are presented in Figure 6. The heating and cooling technologies work with the electricity or gas provided by the supply side or energy storage technologies, and the provided energy is fed to the energy distribution system to be delivered to each zone. Three kinds of energy systems can be modeled based on the technologies mentioned in Section 2, Section 3 and Section 4 to protect the DHW demand, space heating demand and space cooling demand. They can be summarized as follows:

• Electricity system: The onsite energy generation technologies, including PV, WT and CHP, produce electricity for the heating or cooling system and home appliances of the building. The surplus electricity can be stored in the battery or sold to the grid;
• Heating system: The onsite heating production can be provided by renewable generation; for example, STP, or other environmentally friendly technologies, such as micro-CHP and heat pumps. Some traditional heating technologies, such as the electric heater, can also be included in this framework. DHW and space heating are the two main domains of heating demand in the building. Thermal heat can be distributed to each zone by water-based technologies, such as HR systems. The surplus thermal heat can be stored in HWSTs or PCMs;
• Cooling system: Concerning buildings in hot climates, the space cooling demand can be covered by cooling technologies, such as evaporative cooling, chillers and CCHPs. Analogously, the cold energy can be distributed to each zone of the building with water-based techniques, such as chilled beams. The surplus cold energy can also be stored in PCMs.

The aforementioned heating and cooling technologies are all based on water. They are prevalently used in buildings nowadays due to their high efficiency, and can be connected and integrated according to Figure 6. Moreover, even if this work does not cover all kinds of technologies, it is extendable because other water-based technologies can be modeled similarly and added to this framework. With this kind of integration, the modeling and control of buildings can both be carried out in the mathematical programming environment, so the running time can be significantly decreased with the same level of accuracy. With the principle to make the models as simple as possible while keeping the accuracy, this work can be a foundation for further research in various scenarios, from the sub-system level to the building level.

The combinations of these proposed models can be integrated with different control strategies to set up an effective EMS for a wide range of types of buildings. The control variables can be continuous variables, such as energy flows, part load ratios and setpoint temperatures, or discrete variables, such as on/off states [203,204]. Simple EMS with RBC is preferred in simple scenarios where a single objective and a few technologies are considered. It defines the operation thresholds or rule sets for each technology based on hourly electricity prices or other incentives, and provides acceptable results comparable to advanced strategies [166]. However, RBC is not so attractive due to its inferior energy performance in complex cases when multiple and competitive objectives need to be considered and more technologies are installed. MPC formulation is preferred in these complex cases, since a larger saving potential can be achieved by solving an optimization problem [186]. MPC can search for the trade-off between different objectives but is costly to implement. Data-driven methods are sometimes included in the process of setting predictive models and solving optimization algorithms to reduce the calculation burden [189].

Discussing different components of an EMS framework, this paper is valuable to the research of several open problems in this field. One of them could be handling thermal comfort when designing EMS. Thermal comfort explains the human satisfactory sense of the thermal environment. Thermal comfort refers to a number of conditions in which the majority of people feel comfortable. To monitor and control the performance of a building, the thermal comfort should be maintained all the time. Hence, the quantification of it is necessary for control law design and optimization. Covering all of the thermal comfort conditions is outside the scope of this review; however, the most popular ones often considered in EMS are presented. There are two main types of factors that significantly affect thermal comfort, categorized as environmental factors and personal factors, respectively [177,178]. The environmental factors include the indoor air temperature, indoor operative temperature, indoor air relative humidity and indoor air velocity. The personal factors include the clothing level and metabolic heat. Each factor has its influence on thermal comfort, and thus all factors in the two categories must be taken into account.

The Fanger method is one of the most popular way to quantify thermal comfort: it defined the so-called PMV index [176] by using the abovementioned factors, and its formula can be found in [205]. The PMV index is expressed based on a seven-point thermal sensation scale (−3, −2, −1, 0, +1, +2 and 3). EN 15251 [206] describes different categories and different levels for the occupants and the comfort conditions. According to ISO 7730 [207], an ideal PMV value ranges between −0.5 and 0.5, where a null PMV value (PMV = 0) represents the best thermal comfort.

Variations in PMV have a direct influence on the level of the building performance [208]. For example, it is shown that up to 33.5% of energy can be saved by a PMV-based cooling system [209]. Apart from energy-saving, the PMV-based control can also achieve a better thermal comfort compared with the temperature-based control [210]. Meanwhile, PMV can be coupled with both RBC and MPC, leading to an approximately 40% energy reduction, as shown in a commercial building in Florida, USA [211].

Another open problem could be the unified framework for building the EMS design based on the connection of different sections, as shown in Figure 7. This framework can be surely constructed based on the presented mathematical models of assemblable buildings physics, energy systems and energy storage, in connection with thermal comfort and EMS based on different control strategies. Another open problem worth mentioning is the application of artificial intelligence techniques to model the different components in Figure 7. One may notice that recent surveys dedicated to this topic either focus on specific equipment [212,213] or they have a more general scope [214,215] but do not explain how to use mathematical models to decrease the complexity of the calculation with data-driven methods.

To conclude, this work has reviewed the mathematical models of building physics and popular technologies, and has presented assemblable models. The fact that the models can be assembled according to different configurations helps the design of an EMS framework that can be used in different locations and climates, and is expandable for various energy systems and control attitudes. The possible connections of these models with thermal comfort and data-driven aspects put forward a few open problems in the study of improving building energy efficiency.