Simplified Building Thermal Model Development and Parameters Evaluation Using a Stochastic Approach

Abstract

This paper proposes an approach to develop building dynamic thermal models that are of paramount importance for controller application. In this context, controller requires a low-order, computationally efficient, and accurate models to achieve higher performance. An efficient building model is developed by having proper structural knowledge of low-order model and identifying its parameter values. Simplified low-order systems can be developed using thermal network models using thermal resistances and capacitances. In order to determine the low-order model parameter values, a specific approach is proposed using a stochastic particle swarm optimization. This method provides a significant approximation of the parameters when compared to the reference model whilst allowing low-order model to achieve 40% to 50% computational efficiency than the reference one. Additionally, extensive simulations are carried to evaluate the proposed simplified model with solar radiation and identified model parameters. The developed simplified model is afterward validated with real data from a case study building where the achieved results clearly show a high degree of accuracy compared to the actual data.

1. Introduction

In the current situation, buildings in EU countries consume more primary energy than other sectors such as: transportation and industry [1]. Rapid increase in the need of new buildings due to the growth on population growth and mass migration to the urban areas accounts for increased energy demand in the buildings. Making all new buildings as nearly zero-energy buildings (NZEB), and retrofitting them structurally and operationally will have a greater impact on the reduction of primary energy consumption [2,3,4].

There has been growing interest in the development of multi-objective controllers [5,6] to improve the overall performance of buildings. A Building Automation System (BAS), or a Building Automation Control System (BACS), is a intelligent controller aiming to obtain comfort, and energy efficiency in the building [7]. Indoor comforts and energy consumption are affected by various parameters, particularly occupants behavior [8], construction type [9], and weather conditions [10]. The occupants behavior and environmental conditions are volatile in nature. They cause unexpected changes in indoor comfort and could result in excessive energy consumption [11,12]. Hence, an intelligent controller should be able to consider occupancy behavior along with weather conditions to predict and control indoor comfort while optimizing energy consumption. A promising controller for BACS is Model Predictive Control (MPC). MPC applied to Heating, Ventilation, and Air-Conditioning (HVAC) has shown excellent performance abilities to achieve comfort and energy optimization in buildings. MPC has the ability to make a better trade-off between indoor comfort and energy consumption of buildings [13] while considering constraints (occupants behavior and weather conditions).

The performance efficiency of MPC applied for HVAC mainly depends on accuracy and computational efficiency of the building model [14,15,16]. The building model (zone model) should replicate dynamics of the heat transfer in building, able to predict zone heat loads/air temperatures as a function of controllable (HVAC systems) and uncontrollable measurable inputs (weather, internal heat gains, solar gains, etc.), and simplified enough to connect to optimization techniques.

Three fundamental solutions to develop simplified models are data-driven method, software-based model and resistor-capacitor (RC) thermal networks. MPC developed based on thermal RC network models have been successful in controlling the HVAC systems in building [16]. This study focuses on development of building model using thermal RC networks and its parameters identification.

Purpose of developing a thermal RC network method is to realize a reduced order model to represent the building thermal behavior [15]. Here, the idea is to represent a thermal system as a linear electrical circuit using the thermal-electrical analogy principle. Then the thermal network model represented as electrical circuit can be solved using state-space equations. Thermal - Electrical system analogy is given in

Table 1. Thermal to Electrical system analogy.

	Thermal System	Electrical System
Source	Temperature (T)	Voltage (V)
	Heat flux ($\varphi$)	Current (I)
Element	Thermal conductivity (k)	Conductivity ($\sigma$)
	Thermal resistance (R)	Electrical resistance (R)
	Thermal capacity (C)	Electrical capacitance (C)

.

The equivalent thermal network circuit of a zone is determined by incorporating models of building envelope such as: walls, windows, and internal mass, etc. In thermal RC network model, the building is split into a network of nodes with interconnecting paths, through which the energy flows [17,18]. The application of this method differs primarily based on the choice of nodes on which energy balance is applied. These model can be developed for two categories:

• Model for a building envelope (walls, floors, roofs, etc.). These models are then used to develop complete zone model.
• Model for a complete zone.

Whilst, there is not huge difference between two methods, the model for a complete zone is developed by aggregating the individual envelope models into a single model. Whereas, the specific model for building envelope are significant in knowing surface temperatures, analysis of individual envelope dynamics, etc. The building envelop thermal network models are normally represented by 2R1C or 3R2C or 4R2C networks [19]. Furthermore, the parameter (resistors and capacitors) values of the thermal network model has significant impact on accuracy of the model performance.

These parameters identification can be achieved by two ways: analytical [19,20] and numerical approach. In analytical method, parameter values is obtained by solving set of algebraic equations and approximation of parameters’ value by optimization techniques in numerical methods. Advantage of the numerical approaches is that they have higher accuracy as a result of the optimization process. Whereas, analytical approaches have limited applicability for few combinations [21]. Xu and Wang [22], introduced an approach to identify parameter values by using genetic algorithm. Fraisse et al. [19], proposed a technique to evaluate resistors and capacitors values analytically, then analysed the performance of various RC combinations (2R1C, 1R2C, 3R2C, and 3R4C). These parameter values can be used to develop a reduced order model to realize thermal behavior of buildings while considering radiative and convective loads from the internal and external sources. In many previous studies, simplified model approaches focus on developing the system model of reduced system of equations which are computationally faster than the conventional models [23]. However, every change in the definition of the associated energy systems, can lead to a totally different model. Thus, the model simplification procedure has to be performed again in this case.

This paper proposes a method for parametric identification of thermal network model (3R2C) by using particle swarm optimization (PSO) algorithm [24]. A simplified model is then developed these identified values for predicting the thermal dynamics of a zone. A reference model is developed using well established Crank-Nicolson finite difference method to represent the energy flow through the building envelope [25]. In literature, parameters are identified by comparing simplified model with the reference model response for a step excitation of input variables. However, in reality the input variables vary dynamically, for example the output temperature varies periodically. Therefore, in this paper, a step and a periodic input excitation given and responses are then compared between reference model and second-order thermal network models. Subsequently, PSO algorithm is applied to minimize the root mean square error (RMSE) between the reference and thermal network models. The calibrated model with the optimized parameter values has actual dynamics of the building envelope, and 40% to 50% computationally faster than the reference model. Using this approach a simplified model developed for a classroom in case study building (container building) and predicted indoor temperature of the model is validated against the real data obtained from the container building.

2. Parameters Identification

The proposed method involves the reference model development for second-order thermal model validation and its parameters identification. Three types of major heat transfers occur in building envelope are conduction, convection, and radiation (see Figure 1). The effects of these heat transfers to be considered while developing the model. The conduction heat transfer responsible for higher percentage of heat transfer compare to the other two processes.

2.1. Reference Model

The conduction heat transfer in a wall under steady-state conditions can be defined using Fourier’s law of heat conduction [26]:

$${\dot{Q}}_{cond}=-kA\frac{dT}{dx}\left(W\right)$$

(1)

where ${\dot{Q}}_{cond}$—conduction heat transfer through building envelope (W),

k—thermal conductivity of the material (W/m·K),

T—temperature (K),

A—surface area of the wall (m${}^2$), and

x—spatial coordinate (m)

The thermal energy flow through multilayer homogeneous wall is shown in Figure 1. Equation for the heat transfer through building envelope can be obtained by making these assumptions:

• the conduction heat transfer is considered to be one-dimensional due to the high ratio of height to thickness of the building envelope. This lead to maximum heat transfer in one direction,
• there is no heat source or sink within the wall,
• effects of thermal bridge are neglected, and
• properties of the building materials are independent of temperature.

The temperature at any given point is the function of space x and time t, i.e., $T=T(x,t)$. By making above assumptions the heat transfer Equation (1) becomes:

$$\frac{\partial T(x,t)}{\partial t}=\alpha \frac{{\partial }^{2}T(x,t)}{\partial x^2}\mathrm{for}0<x<L,t>0,$$

(2)

where $\alpha =\frac{k}{\rho C_p}$—thermal diffusivity (m${}^2$/s),

L—thickness of the wall (m),

k—thermal conductivity of the material (W/m·K),

$\rho$—density (kg/m${}^3$),

$C_p$—heat capacity of the material (J/kg·K).

Equation (2), requires to be solved with boundary conditions at two surfaces of the envelope, and initial conditions. Due to the effect of movement of air near the wall surfaces results in convective heat transfer, hence convective boundaries are formulated as follows:

$$\begin{array}{c}\hfill k{\left(\frac{\partial T}{\partial x}\right)}_{x=0}=h_c[T_s\left(t\right)-T_{x=0}\left(t\right)]\\ \hfill k{\left(\frac{\partial T}{\partial x}\right)}_{x=L}=h_c[T_{x=L}\left(t\right)-T_i\left(t\right)]\end{array}$$

(3)

The partial differential equation (PDE) (2), along with convective boundary conditions (3), can be solved analytically, by using well established methods: Laplace transforms or Separation of variables or others. However, heat transfer through building envelope is time-dependent, and is generally difficult to model using analytical approach. The finite difference Crank-Nicolson method is therefore implemented to approximate solutions at finite time and space, due to its proven ability to solve PDEs [25].

Crank-Nicolson Finite Difference Model (CNFDM)

In order to formulate heat transfer model using CNFDM , the multilayer wall of Figure 1 is discretized into equal segments of spatial width $\Delta x$ along its thickness. Therefore a set of algebraic equations are developed by discretizing the governing Equations (2) and (3) corresponding boundary conditions, using unconditionally stable CNFDM. The general form of CNFDM is detailed as follows:

$$\begin{array}{cc}\hfill T_i^{(t+\Delta t)}-T_i^{\left(t\right)}=& \mathsf{\Psi }[(T_{i-1}^{(t+\Delta t)}-2T_i^{(t+\Delta t)}+T_{i+1}^{(t+\Delta t)})+(T_{i-1}^{\left(t\right)}-2T_i^{\left(t\right)}+T_{i+1}^{\left(t\right)})]\hfill \end{array}$$

(4)

where $\mathsf{\Psi }$ = $k\Delta t/2\rho C\Delta x^2$,

Equation (4) applied for interior nodes is expressed as:

$$\begin{array}{c}\hfill (2+2{\mathsf{\Psi }}_l)T_i^{(t+\Delta t)}=(2-2{\mathsf{\Psi }}_l)T_i^{\left(t\right)}+{\mathsf{\Psi }}_l(T_{i-1}^{\left(t\right)}+T_{i+1}^{\left(t\right)}+T_{i-1}^{(t+\Delta t)}+T_{i+1}^{(t+\Delta t)})\end{array}$$

(5)

where ${\mathsf{\Psi }}_l$ = $k_l\Delta t/2{\rho }_l{C}_l\Delta x^2$ for sublayers $l=1,2,3,4\dots n$, discretized equal segments of the wall. Nodes at the boundary of two intermediate layers (node between $i=2$ and $i=3$ in Figure 1), the thermo-physical properties of these two material layers have an impact on the transient heat transfer. The corresponding equations from CNFDM approach are as follows:

$$\begin{array}{cc}\hfill (2+{\mathsf{\Psi }}_{l1}+2{\mathsf{\Psi }}_{l2})T_i^{(t+\Delta t)}=& (2-{\mathsf{\Psi }}_{l1}-2{\mathsf{\Psi }}_{l2})T_i^{\left(t\right)}+{\mathsf{\Psi }}_{l1}{T}_{i-1}^{\left(t\right)}+{\mathsf{\Psi }}_{l2}{T}_{i+1}^{\left(t\right)}+{\mathsf{\Psi }}_{l1}{T}_{i-1}^{(t+\Delta t)}+{\mathsf{\Psi }}_{l2}{T}_{i+1}^{(t+\Delta t)}\hfill \end{array}$$

(6)

where ${\mathsf{\Psi }}_{l1}=\frac{{\lambda }_{l1}\Delta t}{({\rho }_{l1}{C}_{l1}+{\rho }_{l2}{C}_{l2})\Delta x^2}$ and ${\mathsf{\Psi }}_{l2}=\frac{{\lambda }_{l2}\Delta t}{({\rho }_{l1}{C}_{l1}+{\rho }_{l2}{C}_{l2})\Delta x^2}$.

The convection boundary conditions at node $x=0$ and $x=L$ are formulated as listed below:

$$\begin{array}{cc}\hfill (2+2\mathsf{\Psi }+2H)T_1^{(t+\Delta t)}=& (2-2\mathsf{\Psi }-2H)T_1^{\left(t\right)}+2\mathsf{\Psi }(T_2^{\left(t\right)}+T_2^{(t+\Delta t)})+H(T_e^{\left(t\right)}+T_e^{(t+\Delta t)})\hfill \end{array}$$

(7)

$$\begin{array}{cc}\hfill (2+2\mathsf{\Psi }+2H)T_N^{(t+\Delta t)}=& (2-2\mathsf{\Psi }-2H)T_N^{\left(t\right)}+2\mathsf{\Psi }(T_{N-1}^{\left(t\right)}+T_{N-1}^{(t+\Delta t)})+H(T_e^{\left(t\right)}+T_e^{(t+\Delta t)})\hfill \end{array}$$

(8)

where $H=h_c\Delta t/\rho C\Delta x$,

$T_e$ = ambient temperature for boundary at $x=0$, and zone air temperature for boundary at $x=L$.

The set of CNFDM equations of energy conservation can be expressed in tridiagonal matrix.

$$A{T}_{1..n}^{(t+\Delta t)}=B{T}_{1..n}^{\left(t\right)}$$

(9)

where A is the matrix with future values and B is the matrix with present values with respect to time coefficients $(t+\Delta t)$ and $\left(t\right)$, respectively. T is the temperature vector at $\left(t\right)$ and $(t+\Delta t)$. By applying Thomas algorithm with the given initial and boundary conditions, the future temperature values are determined [27].

2.2. 3R2C Thermal Network Model

The thermal network model approach is used to simplify the complex modeling of building dynamics, and to obtain reduced order simplified models of the buildings based on state-space equations. The 3R2C network model consists of three R’s and two C’s representing resistances and capacitances of composite homogeneous wall, respectively. Lorenz and Masy [28] proposed among the early implementations of a thermal network model of configuration 2R1C (two resistors and one capacitor). Gouda et al. [21] proposed improvised model of configuration 3R2C to obtain accurate models by non-linear optimization. The paper concluded that 3R2C model with proper parameter values is best suited for practical applications in comparison with 2R1C model.

Generally, thermal network models for envelope are modeled as 2R1C or 2R2C or 3R2C network [29,30,31]. Nevertheless, application of 3R2C model can be seen regularly in the literature due to its performance accuracy and low computational cost. In 3R2C model (see Figure 2), $R_1$, $R_2$, and $R_3$ are resistors of wall thermal resistivity, and $C_1$ and $C_2$ are capacitors representing thermal capacitances of the composite homogeneous wall. The degree of differential equations is equal to the number of capacitors in the network. Consequently, second-order differential Equations (10) and (11) are obtained from 3R2C model.

$$\frac{\mathrm{d}{T}_{c1}}{\mathrm{d}t}=\frac{T_{out}}{R_1{C}_1}-\frac{T_{c1}}{R_1{C}_1}-\frac{T_{c1}}{R_2{C}_1}+\frac{T_{c2}}{R_2{C}_1}$$

(10)

$$\frac{\mathrm{d}{T}_{c2}}{\mathrm{d}t}=\frac{T_{c1}}{R_2{C}_2}-\frac{T_{c2}}{R_2{C}_2}-\frac{T_{c2}}{R_3{C}_2}+\frac{T_{in}}{R_3{C}_2}+\frac{Q_{in}}{C_2}$$

(11)

These set of equations can be expressed using state-space formulation:

$$\begin{array}{cc}\hfill \dot{T}& =AT+BU\hfill \\ \hfill y& =CT+DU\hfill \end{array}$$

(12)

where T is the temperature vector, U is the input vector, y is the output vector, A is the state matrix coefficients related to state vector, B is the input matrix values related to input vector, C is the output matrix values related to state vector, and D is the direct transition matrix.

2.3. Particle Swarm Optimization

Composite multilayer walls have different thermo-physical properties for each layer, thus setting the values of resistors and capacitors of a thermal network model to replicate dynamics of the reference model is a difficult task. As these values have significant impact on performance of the model, it is essential to obtain correct values for each resistor and capacitor of a 3R2C model that emulate the real dynamics of a building envelope.

In this study, a PSO algorithm is applied for parameter values identification of a thermal network model. PSO is a population-based stochastic algorithm [32,33,34], and is suitable for discontinuous non-linear systems with convergence behavior and robustness as key features. The particles of PSO update themselves with the velocity, they even have a memory of the preceding location. In addition, PSO particles have individual best position and move towards global best position by sharing information with other particles. This characteristics makes PSO to converge faster compared to other algorithms for example GA (crossover) [35]. Therefore, PSO can be found well suited to parametric optimization processes when compared to other evolutionary algorithms.

A constrained PSO algorithm was developed to identify model parameter values by minimizing the RMSE between reference and 3R2C thermal network models. The RMSE is minimized by using objective function (13).

$$minf(x1,x2,x4)=\sqrt{\frac{{\sum }_{k=1}^n{(Tfd{m}_k-Tsi{m}_k)}^2}{n}}$$

(13)

$$\begin{array}{cc}\hfill \mathrm{subject}\mathrm{to}\mathrm{constraints}:& \hfill \\ \hfill & x_1+x_2+x_3=1;x_1,x_2,x_3>0\hfill \\ \hfill & x_4+x_5=1;x_4,x_5>0;\hfill \\ \hfill & 1-(x_1+x_2)>0;1-x_4>0\hfill \end{array}$$

(14)

R and C values of 3R2C model are obtained from below expressions:

$$\begin{array}{cc}\hfill & R_1=x_1{R}_{total},R_2=x_2{R}_{total},R_3=x_3{R}_{total},\hfill \\ \hfill & C_1=x_4{C}_{total},C_2=x_5{C}_{total}\hfill \\ \hfill & R_{total}=\sum _{j=1}^M{R}_{layer,m},C_{total}=\sum _{j=1}^M{C}_{layer,m}\hfill \end{array}$$

(15)

where $R_{total}$ is the total resistance and $C_{total}$ is the total capacitance values of composite wall, $f(x_1,x_2,x_4)$ = fitness function to be minimized, $Tfd{m}_k$ and $Tsi{m}_k$ = step and periodic input response temperature at interval of (k).

A number of random swarm particles are generated to initialize the search for an optimal fitness of the PSO algorithm. These particles are updated by their personal best ’$pbest$’ and generation’s global best ’$gbest$’ values at each iteration. The particles values are updated based on the equations:

$$\begin{array}{cc}\hfill V_{i,j}^{k+1}=& w{V}_{i,j}^k+c_1{r}_1(Pbes{t}_{i,j}^k-X_{i,j}^k)+c_2{r}_2(Gbes{t}_j^k-X_{i,j}^k)\hfill \\ \hfill X_{i,j}^{k+1}=& X_{i,j}^k+V_{i,j}^{k+1}\hfill \end{array}$$

(16)

where $Pbes{t}_{i,j}^k$ = personal best jth generation of ith particle,

$Gbes{t}_j^k$ = global best of jth generation,

$V_{i,j}^{k+1}$ = updated velocity of the particle,

$X_{i,j}^{k+1}$ = updated position of particle.

The constrained PSO is influenced by a number of parameters, namely number of particles, acceleration coefficients, inertia weight, and number of iterations. These parameter values were selected based on [36].

2.4. Simulation Results of Parametric Identification

In this study 3 types of building envelope datasets are used to demonstrate the suitability of proposed parametric identification approach. These datasets are light, medium, and heavy weight/thermal mass composite walls. Low thermal mass construction walls typically consist of either steel or timber layers with insulating materials. Typically, low mass wall systems are at least partly pre-assembled off-site this indicates that they are pre-fabricated in a industry/factory. With this option we can construct building quicker than the heavy mass walls, as the construction work is less weather dependent. Low mass walls have lower amount of heat storage ability. Whereas, heavy mass walls generally built of precast concrete or concrete blocks. Properly insulated heavy mass walls can store significant amount of heat from sun, thus minimizing cooling and heating load while enhancing indoor comfort. Many of the conventional buildings are built using heavy mass walls. New buildings are recently using low mass walls. Therefore, the parameters identification approach is applied to different thermal mass walls.

The multilayer wall layer’s thermo-physical properties and wall composition details are provided in

Table 2. Multilayer wall classification, thermal properties, and optimized R and C values.

Construction Class	Thermal Properties				Parametric Values R (m${}^2$·K/W and C (kJ/m${}^2$·K)
	Thickness	Conductivity	Density	Specific Heat
	mm	W/(m·K)	kg/m${}^3$	kJ/(kg·K)	${\mathit{R}}_{\mathit{total}}$	${\mathit{C}}_{\mathit{total}}$	${\mathit{R}}_1$	${\mathit{R}}_2$	${\mathit{R}}_3$	${\mathit{C}}_1$	${\mathit{C}}_2$
Light-weight (LW)					3.1498	76.852	0.29477	2.7812	0.07383	20.694	56.157
Stucco	25.00	0.692	1858	0.84
Insulation (batt)	125.00	0.043	91	0.96
Plaster/Gypsum	20.00	0.727	1602	0.84
Medium-weight (MW)					3.8238	183.724	0.0937	3.6735	0.0565	69.664	114.059
Brick	101.60	0.89	1920	0.79
Insulation board	50.80	0.03	43	1.21
Air space	50.00	0.02514	1.205	1.00
Gypsum	20.00	0.727	1602	0.84
Heavy-weight (HW)					2.1917	402.102	0.1417	1.9018	0.1481	205.196	196.906
Brick	101.60	0.89	1920	0.79
Heavyweight concrete	203.2	0.53	1280	0.84
Insulation board	50.80	0.03	43	1.21
Gypsum	20.00	0.727	1602	0.84

. These thermo-physical properties and composition of layers data are obtained from the ASHRAE Handbook of Fundamentals [37]. In order to identify the parameter values of thermal network model and to validate the its dynamics, a reference model is developed. In addition, two other thermal network models (model (I) and (II)) are developed for comparison purposes. The resistors and capacitors values of model (I) are assigned in such a way that the total resistance and capacitance values are allocated equally for $R_1$, $R_2$, $R_3$, $C_1$, and $C_2$ parameters. Whereas, the values of resistors and capacitors of model (II) are assigned same as the building envelope layer’s thermal resistance and capacitance meaning that the layer (outside) value is assigned to resistor $R_1$. Capacitors $C_1$ and $C_2$ values are allocated in the same way.

Furthermore, the CNFDM-based reference model is divided into 80 equal segments to have higher accuracy with initial conditions $T(x,t)=0$ at $t=0,\forall x\in [0,L]$, boundary conditions $T(x,t)=u\left(t\right)$$\forall t>0,$ at $x=0$ and $x=L$, $T(x,t)=T_e\left(t\right)$, $\forall t>0$. A unit step and periodic input is given to the models (optimized 3R2C model, model (I), and model (II)) and output indoor temperature is measured.

The minimization of the objective function by PSO algorithm is shown in Figure 3. The minimized fitness function value is plotted against the particle generations in which each generation consists of 100 particles. Initial particles are randomly generated within the bounded values, hence few particles in first few iterations have violated the given constraints, those are omitted in the plot.

The developed models were simulated in Python using the following computational facilities: Intel Core i5-7300U, CPU 2.60 GHz and 8 GB (RAM) under operating. Accordingly, the average execution time of PSO algorithm is 30.33 seconds (1000 maximum iterations). The 3R2C thermal network model was around 40% to 50% computationally faster over reference model.

The comparison between temperature responses of heat conduction transfer of the reference model and the other 3 thermal network models for different thermal mass multilayer walls is represented in Figure 4, Figure 5 and Figure 6. For low thermal mass wall, the results in Figure 4 show that the temperature response of optimised model has great fitness to the reference model characteristics. However, the model (II) has acceptable accuracy but the model (I) has taken more time to reach stability with big difference in response compared to the reference model response, hence produces large errors in real applications when there is dynamic change in temperature. Meanwhile, the model (II) of the low thermal mass wall has shown relatively good accuracy with reference model. The parameter values of the model (II) are same as of the 3 layers in the low mass wall. The results show that for low thermal mass model (II) can be applied in real time applications due to its acceptable performance. While, in the case of accuracy is very important, optimized model is the best solution. The root mean square error (RMSE) for different construction class walls are given in

Table 3. Root mean square error (RMSE) between reference and simplified thermal network model.

Construction Class	Step Excitation	Periodic Excitation
Light-weight(LW)	1.239 × ${10}^{-3}$	1.651 × ${10}^{-2}$
Medium-weight(MW)	1.311 × ${10}^{-3}$	1.365 × ${10}^{-2}$
Heavy-weight(HW)	1.495 × ${10}^{-3}$	2.208 × ${10}^{-2}$

.

In the case of medium and heavy thermal mass walls, they have 4 layers, and the parameter values of model (II) are assigned as follows: first layer resistance value to $R_1$, combined resistance value of second and third layer to $R_2$, and third layer value to $R_3$. The temperature response of the model (I) has poor fitness with reference model for these two walls. In Figure 6, the model (II) of heavy thermal mass wall has followed the trajectory of reference model in both step excitation and periodic excitation but has poor fitting with the medium thermal mass wall (see Figure 5). Furthermore, the temperature response of model with optimized parameter values has closely followed the reference model response dynamics in both cases of step and periodic excitation. In this context, the proposed 3R2C model with optimal parameter values is strongly recommended because of its increased accuracy and computational efficiency.

Figure 4, Figure 5 and Figure 6 shows that the time-lag difference between the different thermal mass walls, the time response of the model is less in low mass wall compare to other two walls. The time-lag in the heavy mass wall is almost 2 days to reach the steady state.

Comparison with Conduction Transfer Function (CTF) Model

There are various methods to simulate the performance of a building. These methods can be put into two categories: (1) Analytical, and (2) Numerical. In this paper, numerical finite difference method has been chosen to develop reference model. Many simulation tools are developed based on the finite difference approach [38]. However, many recent well established simulation tools (EnergyPlus, TRNSYS, etc.) use analytical conduct transfer function (CTF) method to perform calculations.

The parameters of simplified thermal network model are identified by comparing it with the finite difference model. Furthermore, a comparison simulation has been conducted to compare thermal network performance with the analytical conduction transfer function method. The CTF method was introduced by Mitalas and Stephenson [39]. A transient conduction calculations performed on multilayer wall using CTF method is detailed in [40], all coefficients for this wall was provided along with cooling load output, and input outdoor and indoor temperatures. Applying the parameters identification technique detailed in Section 2 on example wall, the resistor and capacitor values are obtained. The simplified thermal network model is then developed for the example wall and compared with the CTF model in Figure 7. A very good agreement can be noticed between optimized thermal network and CTF model. The symmetric mean absolute percentage error (sMAPE) between the both models was found to be 0.02008 (≈20%). The developed parametric identification method based on CNFDM, presents very good prediction accuracy when compared with CTF model.

3. Detailed Modeling Approach for a Case Study Building

3.1. Building Description

The CESI smart container building is a multi-purpose building, the rooms are used as classrooms and laboratories, located in the CESI Campus of Nanterre, France. This building is constructed under the program French Programmes d’Investissement d’Avenir (PIA-France). The building meets the French energy standard BBC (Bâtiment à Basse Consommation) [41].

Figure 8a shows the CESI smart building (Nanterre, France). It is an intelligent and connected building built in partnership with companies CISCO, Philips Lighting and Vinci Energies. The 200 m² building is made on the basis of 16 recycled maritime containers. The SMART Building integrates different technologies and data sensors to detect various physical parameters. These include LED lights powered and controlled by POE (Power Over Ethernet), opening sensors, POE temperature and humidity sensors in each room, POE cameras and the presence of a weather station on the roof. Thus, the demonstrator allows the collection of data on brightness, temperature, humidity, presence, sash opening, energy consumption and occupancy. The server is capable of controlling the systems and modifying the thermal and luminous ambience’s. The heating and cooling energy is supplied to indoors through installed heat pump and ventilation system. The building is south facing with a door opening, and south wall is almost covered with windows. A classroom was selected as a case study room (see Figure 8b). The objective is to develop a simplified thermal network model using the parametric identification and validate it against the measured temperature values. Similar approach can be applied to other rooms as they are similar to the case study one.

3.2. Thermal Network Model—CESI Smart Building

The heterogeneous nature of buildings and parameters that influence its performance makes modeling of a zone (case study room) highly complex. Particularly, modeling of building using analytical approach is infeasible because of its non-linear behavior.

The thermal network approach discretize the complex building into multiple zone, where each zone is assumed to have properly mixed well and the building walls have uniform temperatures throughout its volume. The zones are then modeled by means of the network of nodes with interconnecting heat flow paths. The heat transfer can occur through conduction, convection, and radiation. Heat gains from different means: internal, solar radiation, etc,. are lumped in the thermal network nodes, and the heat storage in building thermal mass is in thermal capacitances. Temperature and heat flows between the nodes are determined by energy balance approach as described in Section 2.2. The resultant formulations are a set of coupled ordinary differential and algebraic equations that can be solved using the state space equation.

The simplified thermal network model was developed by making the following assumptions:

• heat conduction occurs through the building envelope,
• convective heat transfer at building envelope surfaces and floors,
• solar gains through windows and solar radiation absorption in external walls,
• radiation heat transfer within the zone walls,
• effect of thermal bridge is neglected,
• heat storage is not considered in windows, and
• impact of wind velocity variation on the convective heat exchange coefficient of the wall surface is neglected, hence the convective resistances are considered constant.

The CESI container building external walls are built using variety of materials for different sides: east, west, north and south. Therefore, 3R2C model is developed individually for each side of external walls, floor, and roof.

Considering all the above hypothesis, a thermal model is developed with 6 sets of 3R2C networks and the equations obtained by energy balance at each node from the schematic shown in Figure 9:

$$\begin{array}{cc}\hfill C_z\frac{d{T}_z}{dt}=& \sum \frac{T_{1j}-T_z}{R_{1j}}+\frac{(T_{amb}-T_z)}{R_{win}}+{\dot{Q}}_{in}\hfill \end{array}$$

$$\begin{array}{cc}\hfill & +{\dot{Q}}_{s,wi}+{\dot{Q}}_{heat}\hfill \end{array}$$

(17)

$$\begin{array}{cc}\hfill C_{1j}\frac{d{T}_{1j}}{dt}=& \frac{T_{2j}-T_{1j}}{R_{1j2}}+\frac{T_{in-1j}-T_{1j}}{R_{1j3}}\hfill \end{array}$$

(18)

$$\begin{array}{cc}\hfill C_{2j}\frac{d{T}_{2j}}{dt}=& \frac{T_{ou-1j}-T_{2j}}{R_{1j1}}+\frac{T_{1j}-T_{2j}}{R_{1j2}}+{\dot{Q}}_{s,wj}\hfill \end{array}$$

(19)

$$\begin{array}{cc}\hfill {\dot{Q}}_{s,w}=& {\alpha }_w{A}_w{q}_{rad}^{"}\hfill \end{array}$$

(20)

where, $C_z$ is the zone air capacitance, $C_{1j}$ is the wall capacitance indoor side, $C_{2j}$ is capacitance at outdoor side of the wall, $T_{amb}$ is ambient temperature, ${\dot{Q}}_{heat,pump}$ is heating or cooling energy supplied from heat pump, ${\dot{Q}}_{vent}$ is heating or cooling energy supplied from ventilation system, ${\dot{Q}}_{in}$ is heat gains from occupants and electrical appliances, ${\dot{Q}}_{s,wi}$ is solar gains through windows, ${\dot{Q}}_{s,wj}$ is solar gains from the building external walls, and ${\dot{Q}}_{heat}$ is total supplied energy from ${\dot{Q}}_{heat,pump}$ and ${\dot{Q}}_{vent}$.

The first term in right hand side of (17) represents heat transfer through the envelope (external walls), these external walls are exposed to outside temperature and solar radiation. The second term represents heat transfer through the windows, which allow solar radiations directly to the zone. The internal heat gain ${\dot{Q}}_{in}$ is the sum of all internal gains (inhabitants and electrical equipments) and ${\dot{Q}}_{heat}$ is the sum of energy supplied from heating equipments (heat pump and ventilation system). The model output is zone temperature $T_z$, and the influential parameters of zone temperature are: ambient air temperature, internal free gains, solar radiation. Internal gains from the occupants can be exploited to heat the zone. Since, the chosen room is a class room, the students (adults) will be mostly sitting and reading. The assumed heat gains from one adult student is equal to 116 W in summer and 114 W in winter [37,42]. The measurable parameters and data collected from the building include outdoor air temperature, indoor air temperature, global horizontal irradiance, relative humidity, ventilation air flow rate, and occupants attendance data. The measurement ranges and associated error ranges of the installed sensors are listed in

Table 4. Measurement and associated error ranges of sensors.

Sensors	Measurement Range	Error Range
Outdoor temperature	−40 °C to +75 °C	±1 °C
Indoor temperature	−20 °C to +70 °C	±1 °C
Indoor illuminance	50 lx to 20,000 lx	±20%
Relative humidity	5% to 95%	±5%

.

The model consists of overall 13 nodes: two nodes for each envelope wall, floor, roof and a zone node. The installed lighting appliances in the room is low energy consumption hi-tech light systems, thus the gains from these light systems are neglected. However, heat gains from student’s laptop is considered during the simulation. The assumed heat gains from one laptop equalled 30 W [37], and there are no other electrical systems installed in the classroom. Furthermore, the selected controllable input and measurable disturbances are: outside air temperature, internal free gains and solar gains. Heating system is the only controllable input in the model. This thermal network model can be used for real-time HVAC control applications. For example, model predictive controller is applied for comfort and energy optimization based on the simplified thermal network model [13].

The formulation of the thermal network model of the building in state-space representation is given as follows.

$$\begin{array}{cc}\hfill X(k+1)& =AX\left(k\right)+Bu\left(k\right)+B_{1}d\left(k\right)\hfill \\ \hfill y\left(k\right)& =CX\left(k\right)+Du\left(k\right)+D_{1}d\left(k\right)\hfill \end{array}$$

(21)

The model parameters are determined from the approach given in Section 2. In general, some of the measurable but uncontrollable inputs are easily measurable (directly) or readily available concerning to the building location such as: ambient temperature. However, other inputs such as solar irradiation are not readily available or measured as global horizontal irradiance. Hence, a detailed analysis on the solar irradiation incident on the building envelope is essential. To calculate the total irradiation energy on the building envelope, it is needed to obtain radiation energy on each surface.

In the literature, many studies were conducted to calculate the solar irradiation incident on the tilted surface by taking account of the isotropic diffuse sky model. In isotropic diffuse sky model, three components of radiations are considered: isotropic diffuse, diffused radiation reflected from the ground, and beam radiation [43]. The isotropic model is straightforward and calculation of radiation on tilted surfaces becomes simple. However, it fails to take account of circumsolar diffuse and/or horizon brightening components on a inclined surface, this results in the underestimation of solar radiation (around 10 to 15% ). More advanced approaches like Hay–Davies–Klucher–Reindl (HDKR) uses anisotropic sky models, that consist three components: isotropic, circumsolar diffuse and horizon brightening [44,45,46]. Figure 10 shows the distribution of all parts of solar radiation on a tilted surface.

Considering an anisotropic sky model, the incident solar radiation on a tilted surface based on HDKR model, is calculated by [45]:

$$\begin{array}{cc}\hfill I_{Total}=& (I_b+I_d{A}_i)R_b+I_d(1-A_i)\left(\frac{1+cos\beta }{2}\right)\times \left[1+fsi{n}^3\left(\frac{\beta }{2}\right)\right]+I_{gh}{\rho }_g\left(\frac{1-cos\beta }{2}\right)\hfill \end{array}$$

(22)

where $I_b$ = beam radiation,

$I_d$ = diffuse radiation,

$R_b$ = geometric factor,

$\theta$ = angle of incidence,

$\gamma$ = surface azimuth angle,

$\beta$ = angle between tilted surface and the horizontal plane,

$I_{gh}$ = total radiation on horizontal surface,

${\rho }_g$ = ground albedo (0.3).

Anisotropy index ($A_i$) in Equation (22), is a function of the transmittance of the atmosphere for beam radiation. It is the ratio of beam radiation on a horizontal ground surface to extraterrestrial radiation ($I_o$).

$$\begin{array}{c}\hfill A_i=\frac{I_b}{I_o}\end{array}$$

(23)

The modulating factor f in the correction factor is to account for cloudiness.

$$\begin{array}{c}\hfill f=\sqrt{\frac{I_b}{I}}\end{array}$$

(24)

The angles under consideration and other solar angles are illustrated by Figure 10b.

The angle $\theta$ is the incidence angle of the beam radiation on tilted surface, and the angle ${\theta }_T$ is the incidence angle of the beam radiation on the inclined surface. These angles are calculated by:

$$\begin{array}{cc}\hfill R_b=& \frac{cos\theta }{cos{\theta }_z}\hfill \\ \hfill cos\theta =& sin\delta sin\varphi cos\beta -sin\delta cos\varphi sin\beta cos\gamma +cos\delta cos\varphi cos\beta cos\omega \hfill \end{array}$$

(25)

$$\begin{array}{cc}\hfill & +cos\delta sin\varphi sin\beta cos\gamma cos\omega +cos\delta sin\beta sin\gamma sin\omega \hfill \end{array}$$

(26)

$$\begin{array}{cc}cos{\theta }_z=\hfill & sin\varphi sin\delta +cos\delta cos\varphi cos\omega \hfill \end{array}$$

(27)

Since, building envelope walls are perpendicular to the horizontal plane, the $\beta$ angle between the plane of the surface and the horizontal is 90°. The (27), then becomes:

$$\begin{array}{cc}\hfill cos\theta =& -sin\delta cos\varphi cos\gamma +cos\delta sin\varphi cos\gamma cos\omega +cos\delta sin\gamma sin\omega \hfill \end{array}$$

(28)

the total solar radiation on a surface is given by:

$$\begin{array}{cc}\hfill I_{Total}=& (I_b+I_d{A}_i)R_b+I_d(1-A_i)(1+0.3535f)/2+I_{gh}{\rho }_g/2\hfill \end{array}$$

(29)

The total solar radiation for complete building envelope is calculated as follows:

$$\begin{array}{c}\hfill {\dot{Q}}_{sol}=\sum _{j=1}^4{I}_{Ti}{S}_i{\alpha }_i\end{array}$$

(30)

The sums of hourly solar radiation incident on a south facing wall for 5 days in May in Nanterre, France are shown in Figure 11. The solar radiation calculated for isotropic and anisotropic sky models. These two calculated sums of the radiations are compared with the data procured from the weather station on horizontal surface. The ground albedo is considered as 0.2.

3.3. Model Validation

In this section, the developed model will be validated against the measured data of the case study room. Initially, model parameters are identified using the approach given in Section 2. Ambient and room temperatures are measured using the sensors installed in the CESI building, and the solar radiation energy is calculated from the above formulations. The heating input is measured from the heating system installed in the building.

However, the developed model is non-linear in nature because of the radiation heat transfer between the walls, that has fourth power of temperature.

$$\begin{array}{cc}\hfill {\dot{Q}}_{rad}& =ϵ\sigma (T_h^4-T_c^4)\hfill \end{array}$$

(31)

where ${\dot{Q}}_{rad}$ is the radiation heat energy, $ϵ$ is the emissivity coefficient of the wall surface, $\sigma$ is the Stefan-Boltzmann constant, and T is surface temperature of the wall (°K).

Typically, the effect of heat transfer by radiation is either neglected or considered as constant. Since, the temperature variation between two walls is minimal on an absolute scale, i.e., $T\left[K\right]=273+T_{surface}$. Therefore, the effect of change radiative heat energy is almost negligible. Hence, to avoid non-linearity in the model, a constant value is assumed for radiative heat transfer coefficient. This value is determined by linearizing radiative heat transfer near equilibrium using Taylor’s series expansion. This approach of linearizing the model near equilibrium does not generate significant errors, thus the effect on the model is in acceptable range of accuracy [47].

The simplified thermal network model of the zone is validated against the measured data from case study classroom of the CESI LINEACT building. The model is validated for two seasons, two months in summer and winter, respectively.

The developed model is simulated for two months in different seasons, winter (see Figure 12) and summer (see Figure 13). Several key observations are noted here, there is higher variation in the model response for winter compared to the response of summer. This variation can be due to the fact that unpredicted indoor activities, door/windows openings, varying internal gains and air infiltration. In the selected summer period, the building was mostly unoccupied due to the summer holidays. Therefore, there is less variation in the response and in both the model outputs, zone temperature follows the dynamics of the actual measured temperature of the zone. Furthermore, in summer (see Figure 13) the raise in indoor temperature reaching almost 30 °C is due to non-operation of the HVAC system during holidays. The HVAC systems are scheduled to operate only during the office hours. Measured inputs of the model: solar irradiation, occupation, and ambient and adjacent rooms temperatures are shown in Figure 14 and Figure 15. Nevertheless, in both seasons the model reproduces the dynamics of zone temperature of the building efficiently. The results show that the simplified thermal network model with parameter identification is well suited approach to model the building thermal dynamics, and the models can be used for the controller purpose.

4. Conclusions

This paper has proposed an approach to develop simplified thermal network model and its parameters identification. The model parameter values have significant impact its performance. A stochastic PSO algorithm has been adopted to determine the parameters set that provides the optimal approximation of the proposed low-order model dynamics with respect to the reference model dynamics. The constraints for the PSO algorithm were chosen such that the parameter values are not more than the total material resistance and capacitance. The performance of the simplified thermal model with optimized parameters has computational efficiency 40% to 50% higher compared to reference model.

The study results were divided in two sections; firstly the results of simplified model parameters identification presented with the validation against reference model and secondly a complete zone model is developed and validated with the data of case study building. The model validation result has shown a second order thermal network model produce exact thermal behavior with a better accuracy. The developed 3R2C network model is well suited for model-based controller application. Indeed, it has reduced system states number, it is computationally efficient, and can be adapted to any type of buildings.

The modeling approach proposed in this paper will be further extended to model multi-zone buildings, considering inter-zonal heat transfer. Furthermore, it will be used to predict other comfort parameters such as humidity.