# Dynamic models for energy control of smart homes

Christian Ghiaus

► To cite this version:

Christian Ghiaus. Dynamic models for energy control of smart homes. Stephane Ploix; Manar Amayri; Nizar Bouguila. Towards Energy Smart Homes. Algorithms, Technologies, and Applications, Springer, pp.163-198, 2021, 978-3-030-76476-0. 10.1007/978-3-030-76477-7_5. hal-03578578

HAL Id: hal-03578578

https://hal.science/hal-03578578

Submitted on 17 Feb 2022

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Abstract

The dynamic thermal behaviour of the buildings can be represented by thermal networks. These
oriented weighted graphs are graphical representations of the set of equations of heat transfer
by conduction, convection, radiation and advection. Techniques for linearizing the models,
especially the radiative heat exchange, are shown. Being linear (or linearized) models, the
thermal networks can be described by a set of matrices and vectors. An original method for
assembling the thermal networks by merging specified nodes is introduced. The advantage of
assembling over coupling is that the obtained model represents a single set of equations which
can be analysed. Then, it is shown how to transform thermal networks into state-space
representation.

1 Introduction

Dynamic models are widely used for optimization of energy consumption in buildings, both in
design and operation phase. A vast literature exists on modelling and software implementation
for thermal simulation of buildings based of physical laws of energy and mass transfer (Clarke,
2001; EnergyPlus, 2015; TRNSYS 17, 2009).

The design of control systems requires thermal models of the building. Linear Time Invariant (LTI) models, such as state-space representations, transfer functions or zero-pole-gain models, are widely used for this aim. However, obtaining LTI models for the thermal behaviour of the buildings is difficult for two reasons. First, buildings are complex objects modelled by a set of partial differential equations and nonlinear functions that describe the heat transfer by conduction, convection and radiation and energy transport by advection. For the range of variables involved in heat transfer occurring in buildings, these equations can be considered linear or are locally linearizable. The problem is to obtain the models for a whole building by using the models of the components. Second, the models used in heat transfer are thermal networks, while the models used in control theory are LTI models. The problem is to transform models with hundreds of parameters from thermal networks to LTI representation.

This chapter deals with these two problems. It presents two algorithms, one for assembling the
thermal circuits and another for extracting the state-space representation from the thermal
circuit. Although the examples are given for thermal models of buildings, the two algorithms are
of general interest.

2 Thermal networks

There are three forms of heat transfer: conduction, convection and radiation. Heat is also carried by mass flow or advection. Heat transfer can be modelled by thermal circuits or networks composed of nodes and branches. The nodes, which may have capacities, are connected by conductances. There are two type of sources in the thermal networks: temperature (on the branches) and heat flow rate sources (connected to nodes).

<figure><img src="image_2.png" /><figcaption>Figure 1 Sources of: a) temperature; b) heat flow rate.</figcaption></figure>

2.1 Heat sources

2.1.1 Temperature sources

A temperature source (Figure 1a) represents a difference of temperature on a branch which does not change with the flow rate crossing the branch. In the case of buildings, typically they model:

- outdoor air which does not change its temperature with the temperature of the surfaces of the buildings;

- temperature of the ground at a depth at which it is not influenced by the building;

- adjacent spaces that have their temperature controlled by thermostats.

2.1.2 Heat flow rate sources

A heat flow rate source (Figure 1b) represents a thermal energy rate that does not change with the temperature node in which it is entering. In the case of buildings, typically they model:
- solar (or short wave) radiation on the surfaces of the building;

- electrical power from household appliances which is transferred to the indoor air through convection;

- heat transferred through convection from occupants to the indoor air.

The solar (or short wave) radiation in a building is entering through the windows, touches some surface of the indoor walls and then it is multi-reflected. This heat flow rate is modelled by heat flow rate sources. For finding the values of these sources, let us consider radiation exchange between opaque, diffuse, grey surfaces in an enclosure formed by $j = 1 ... n$ surfaces $S_j$ with reflection coefficient $\rho_j$. The view factors between the surfaces are $F_{ij}$. The direct irradiance of each surface is $E_j^0$. The total thermal flux received by the surface $S_i$ directly and after reflection is

$$S_i E_i = S_i E_i^0 + \sum_{j=1}^{n} F_{ji} S_j \rho_j E_j \quad (1)$$

By using the reciprocity relation for view factors, $F_{ij}S_i = F_{ji}S_j$, equation (1) becomes:

$$E_i - \sum_{j=1}^{n} F_{ij} \rho_j E_j = E_i^0 \quad (2)$$

The set of equations (2) can be written in matrix form

$$\begin{bmatrix} 1 - \rho_1 F_{11} & -\rho_2 F_{12} & \dots & -\rho_n F_{1n} \\ -\rho_1 F_{21} & 1 - \rho_2 F_{21} & \dots & -\rho_n F_{2n} \\ \dots & \dots & \dots & \dots \\ -\rho_1 F_{n1} & -\rho_2 F_{n2} & \dots & 1 - \rho_n F_{nn} \end{bmatrix} \begin{bmatrix} E_1 \\ E_2 \\ \dots \\ E_n \end{bmatrix} = \begin{bmatrix} E_1^0 \\ E_2^0 \\ \dots \\ E_n^0 \end{bmatrix} \quad (3)$$

or

$$(I - F\rho)e = e^0 \quad (4)$$

where

$$I = \begin{bmatrix} 1 & 0 & \dots & 0 \\ 0 & 1 & \dots & 0 \\ \dots & \dots & \dots & \dots \\ 0 & 0 & \dots & 1 \end{bmatrix}; F = \begin{bmatrix} F_{11} & F_{12} & \dots & F_{1n} \\ F_{21} & F_{21} & \dots & F_{2n} \\ \dots & \dots & \dots & \dots \\ F_{n1} & F_{2n} & \dots & F_{nn} \end{bmatrix}; \rho = \begin{bmatrix} \rho_1 & 0 & \dots & 0 \\ 0 & \rho_2 & \dots & 0 \\ \dots & \dots & \dots & \dots \\ 0 & 0 & \dots & \rho_3 \end{bmatrix} \quad (5)$$

and the vectors
$$ \mathbf{E} = \begin{bmatrix} E_1 \\ E_2 \\ \vdots \\ E_n \end{bmatrix}; \mathbf{E}^0 = \begin{bmatrix} E_1^0 \\ E_2^0 \\ \vdots \\ E_n^0 \end{bmatrix} \qquad (6) $$

The vector of irradiance of the surfaces is then

$$ \mathbf{E} = (\mathbf{I} - \mathbf{F}\rho)^{-1} \mathbf{E}^0 \qquad (7) $$

Then, the heat flow rate source for each surface $S_i$ is $q_i = E_i S_i$.

Calculating the view factors $F_{ij}$ may be complicated (Bergman, et al., 2011). A simple but rough estimation for rooms with planar surfaces is

$$ F_{ij} \approx \frac{S_j}{\sum_{k \neq j} S_k}; F_{ii} = 0 \qquad (8) $$

which complies with the summation rule of view factors, $\sum_{j=1}^{n} F_{ij} = 1$, but not with the reciprocity relation, $F_{ij}S_i = F_{ji}S_j$. Another simplified estimation,

$$ F_{ij} \approx \frac{S_j}{\sum_k S_k} \qquad (9) $$

complies with the summation rule and with reciprocity relation but results in $F_{ii} \neq 0$.

## 2.2 Heat resistances (or conductances)

The three modes of heat transfer (conduction, convection and radiation) and the heat advection can be modelled by thermal resistances or conductances (Figure 1).
<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /></figure>

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /><figcaption>Figure 2 Conductances in heat transfer and energy advection: a) conduction; b) convection; c) radiation; d) advection.</figcaption></figure>

### 2.2.1 Conduction

Thermal conduction is the heat diffusion in solids in the direction of the temperature gradient (Figure 1a). Fourier law, the equation relating the thermal heat flow rate in a direction $x$, $q_x$, to the temperature gradient, $d\theta/dx$, in the direction $x$ is:

$$q_x = -\kappa S \frac{d\theta}{dx} \qquad (10)$$

where $S$ is the area of the surface perpendicular to the heat flow rate $q_x$. The minus sign shows that heat transfer is from high to low temperature.

<figure><img src="image_7.png" /><figcaption>Figure 3 Steady state thermal conduction: a) stream tube; b) thermal network model.</figcaption></figure>

Let us consider stationary conduction in a stream tube in a homogeneous and isotropic material without internal heat sources (Figure 3a). Since the heat flow rate $q$ is conserved, Fourier law in section s of the streamline is
$$q = -\kappa S \frac{d\theta}{ds} \qquad (11)$$

where the conductivity $\kappa = \kappa(s)$ and the area surface $S = S(s)$ depend on the curvilinear coordinate $s$. By separating the variables, equation (11) becomes:

$$q \frac{ds}{\kappa S} = -d\theta \qquad (12)$$

By integrating the equation (12) from $s_0$ to $s$, the temperature variation with the distance,

$$q \int_{s_0}^{s_1} \frac{ds}{\kappa S} = - \int_{\theta_0}^{\theta_1} d\theta \qquad (13)$$

can be written as

$$qR = \theta_0 - \theta \qquad (14)$$

where

$$R = \int_{s_0}^{s} \frac{ds}{\kappa S} \qquad (15)$$

is the thermal resistance of the stream tube between $s_0$ and $s$.

If there are internal sources, the variation of the heat flow rate along the curvilinear coordinate $ds$ is

$$dq = p \, dV \qquad (16)$$

where $dV$ is the infinitesimal volume. If $ds \to 0$, then $dV = S \, ds + dS \, ds$. By integrating equation (16) between $s_0$ and $s$, $\int_{q_0}^{q} dq = \int_{s_0}^{s} p \, S \, ds$, it becomes:

$$q(s) = \int_{s_0}^{s} p \, S \, ds + q_0 \qquad (17)$$

The flow rate getting out through the surface $S_1$ is:

$$q_1 = \int_{s_0}^{s_1} p \, S \, ds + q_0 \qquad (18)$$

where $q_0$ is the heat flow rate entering through the surface $S_0$. Substituting (18) in (11), we obtain after integration
$$ \theta_1 = \int_{s_0}^{s_1} \frac{1}{\kappa S} \left( \int_{s_0}^{s} -p S \, ds' \right) ds - q_0 \int_{s_0}^{s_1} \frac{1}{\kappa S} ds + \theta_0 \quad (19) $$

By substituting in (19) the expression of thermal resistances given by equation (15),

$$ R_0 = \int_{s_0}^{s_1} \frac{ds}{\kappa S} \quad (20) $$

we obtain

$$ \theta_1 = - \int_{s_0}^{s_1} \frac{1}{\kappa S} \left( \int_{s_0}^{s} p S \, ds' \right) ds - R_0 q_0 + \theta_0 \quad (21) $$

Equations (18) and (21) can be represented by the thermal circuit presented in Figure 3b where the heat rate low source is

$$ f_1 = \int_{s_0}^{s} p S \, ds \quad (22) $$

and the temperature source is

$$ b_0 = \int_{s_0}^{s_1} \frac{1}{\kappa S} \left( \int_{s_0}^{s} p S \, ds' \right) ds \quad (23) $$

Equation (23) can be integrated by parts. By noting $u \equiv \int_{s_0}^{s_1} pSds$ and $v' \equiv \frac{1}{\kappa S}$, the integration by parts $\int uv'ds = uv - \int u'vds$ of equation (23) becomes

$$ b_0 = - \left( \int_{s_0}^{s_1} pSds \right) \left( \int_{s_0}^{s_1} \frac{1}{\kappa S} ds \right) + \int_{s_0}^{s_1} pS \left( \int_{s_0}^{s} \frac{1}{\kappa S} ds' \right) ds \quad (24) $$

By substituting in (24) the expressions of R given by equation (15) and $R_0$ given by equation (20), we obtain:

$$ b_0 = -R_0 \int_{s_0}^{s_1} p S \, ds + \int_{s_0}^{s_1} R p S \, ds \quad (25) $$

With these notations, equation (18) becomes

$$ q_1 = q_0 + f_1 \quad (26) $$

and equation (21) becomes
$$ \theta_0 - \theta_1 + b_0 = R_0 q_0 \qquad (27) $$

where

$$ e_0 = \theta_0 - \theta_1 + b_0 \qquad (28) $$

is the difference of temperature across the thermal resistance $R_0$.

### 2.2.2 Convection

Convection implies heat transfer in fluids (Figure 1b). If steady state is considered, Newton's law of convection is used, which is a phenomenological simplification having an expression similar to Fourier law:

$$ q = hS(\theta_0 - \theta_1) \qquad (29) $$

where $h$ is the convective coefficient determined experimentally (Bergman, et al., 2011). Typical values of the heat convection coefficient are $h_i = 8$ W/m²K inside and $h_o = 25$ W/m²K outside the building (Clarke, 2001).

### 2.2.3 Long wave radiation

Long wave radiation exchange is between two surfaces that are facing each other and that have different temperatures (Figure 1c). The radiative heat flow rate between two black body surfaces i and j is (Bergman, et al., 2011)

$$ q_{ij} = S_i F_{ij} (M_i^o - M_j^o) \qquad (30) $$

where

$S_i$ is the area of surface;

$F_{ij}$ - view factor between surface i and surface j;

$M_i^o$ and $M_j^o$ - the black body radiant emittance of the surfaces i and j, respectively.

By using Stefan - Boltzmann law,

$$ M = \sigma T^4 \qquad (31) $$

where $\sigma$ is Stefan - Boltzmann constant and $T$ is the temperature of the surface expressed in kelvin, equation (30) becomes

$$ q_{ij} = S_i F_{ij} \sigma (T_i^4 - T_j^4) \qquad (32) $$

The two temperatures $T_i^4$ and $T_j^4$ may be linearized around a mean value $\bar{T}$,

$$ T_i^4 = \bar{T}^4 + 4\bar{T}^3(T_i - \bar{T}) \qquad (33) $$
and

$$T_j^4 = \bar{T}^4 + 4\bar{T}^3(T_j - \bar{T}) \quad (34)$$

By subtracting equation (34) from (33), we obtain

$$T_i^4 - T_j^4 = 4\bar{T}^3(T_i - T_j) \quad (35)$$

The exact value of the mean value $\bar{T}$ can be obtained from the equivalence of equation (35) with

$$T_i^4 - T_j^4 = (T_i^2 + T_j^2)(T_i + T_j)(T_i - T_j) \quad (36)$$

as

$$\bar{T} = \sqrt[3]{\frac{1}{4}(T_i^2 + T_j^2)(T_i + T_j)} \quad (37)$$

Substituting (37) in (32), we obtain the linear expression of the heat flow rate:

$$q_{ij} = S_i F_{ij} \sigma 4 \bar{T}^3 (T_i - T_j) \quad (38)$$

where $S_i F_{ij} \sigma 4 \bar{T}^3$ may be considered as a thermal conductance. For $15 \,^\circ\text{C} < \theta < 30 \,^\circ\text{C}$, i.e. $288.15 \, \text{K} < \bar{T} < 303.15 \, \text{K}$, the value of $\sigma 4 \bar{T}^3$ is about 5 W/m$^2$K, more exactly $5.41 < \sigma 4 \bar{T}^3 < 6.31 \, \text{W/m}^2\text{K}$.

The radiosity of a surface represents the radiative fluxes leaving an opaque, diffuse, grey body surface

$$J_i = \varepsilon M_i^o + \rho_i E_i \quad (39)$$

where

$\varepsilon_i M_i^o$ is the emitted radiant flux, with $\varepsilon_i$ the emissivity and $M_i^o$ the black body radiant emittance of the surface $i$;

$\rho_i E_i$ – the reflected radiant flux, with $\rho_i$ the reflectivity and $E_i$ the incident radiant flux on the surface.

The radiative exchange between two opaque, diffuse, grey surfaces may be expressed by a equation similar to (30):

$$q_{ij} = S_i F_{ij} (J_i - J_j) \quad (40)$$

Following the same reasoning as before, the linear expression of (40) is
$$q_{ij} = S_i F_{ij} 4 \sigma \bar{T}^3 (J_i - J_j) \qquad (41)$$

where $\bar{T}$ is given by (37).

<figure><img src="image_8.png" /><figcaption>Figure 4 Radiation networks and their transformation to thermal networks: a) an example of radiation network; b) same network as in (a) but arranged to have a source of emittance on branch; c) thermal network of the radiation network from (b) – the sources and the conductances are changed.</figcaption></figure>

An example of transforming a radiative network into a thermal network is given in Figure 4. The conductances for the radiation network are:

$$G_1 = \frac{\epsilon_1}{1-\epsilon_1} S_1; \quad G_2 = F_{12} S_1; \quad G_3 = \frac{\epsilon_2}{1-\epsilon_2} S_2; \quad G_4 = F_{13} S_1; \quad G_5 = F_{23} S_2 \qquad (42)$$

The oriented incidence matrix **A**, the conductance matrix **G**, and the vector of inputs **b** are:

$$\mathbf{A} = \begin{bmatrix} 1 & 0 \\ -1 & 1 \\ 0 & 1 \\ -1 & 0 \\ 0 & -1 \end{bmatrix}; \quad \mathbf{G} = \begin{bmatrix} G_1 & 0 & 0 & 0 & 0 \\ 0 & G_2 & 0 & 0 & 0 \\ 0 & 0 & G_3 & 0 & 0 \\ 0 & 0 & 0 & G_4 & 0 \\ 0 & 0 & 0 & 0 & G_4 \end{bmatrix}; \quad \mathbf{b} = \begin{bmatrix} M_1^o \\ 0 \\ M_2^o \\ -M_3^o \\ -M_3^o \end{bmatrix} \qquad (43)$$

respectively (see section 3.2 for their definition). The radiosities $j = [J_1 \ J_2]^T$ are given by:

$$\mathbf{j} = (\mathbf{A}^T \mathbf{G} \mathbf{A})^{-1} \mathbf{A}^T \mathbf{G} \mathbf{b} \qquad (44)$$

(see 4.1 for details).

The radiative network, in which the unknowns are radiosities, can be transformed into a thermal network, in which the unknowns in the nodes are temperatures, by using the transformations for emittances and radiosities

$$M^o = \sigma T^4 \quad \text{and} \quad J = \sigma \theta^4 \qquad (45)$$

where $T$ is a temperature source and $\theta$ is an unknown temperature. The conductances for the temperature network become

$$G_1 = \frac{\epsilon_1}{1-\epsilon_1} S_1 \sigma 4 \bar{T}^3; \quad G_2 = F_{12} S_1 \sigma 4 \bar{T}^3; \quad G_3 = \frac{\epsilon_2}{1-\epsilon_2} S_2 \sigma 4 \bar{T}^3; \quad G_4 = F_{13} S_1 \sigma 4 \bar{T}^3; \quad G_5 = F_{23} S_2 \sigma 4 \bar{T}^3 \qquad (46)$$
### 2.2.4 Advection

Energy advection is the transport of energy by a mass flow rate entering and leaving a control volume (Figure 1c). For example, this is the case of energy transported by air in ventilation. The heat flow rate transported is:

$$q_{12} = \dot{m}c(\theta_1 - \theta_2) \qquad (47)$$

where $\dot{m}$ is the mass flow rate and $c$ is the heat capacity of the fluid.

<table><caption>Table 1 Values of airflow rates in air changes per hour as a function of the position of the window (Recknagel, et al., 2007)table 1.12.1-4</caption><thead><tr><td>Position of the window</td><td>Air changes per hour</td></tr></thead><tbody><tr><td>Closed windows, closed doors</td><td>0 to 0.5</td></tr><tr><td>Tilt window, closed jealousy</td><td>0.3 to 1.5</td></tr><tr><td>Tilt window, no jealousy</td><td>0.8 to 4.0</td></tr><tr><td>Half-opened window</td><td>5 to 10</td></tr><tr><td>Full-opened window</td><td>9 to 15</td></tr><tr><td>Windows and French-windows fully opened</td><td>about 40</td></tr></tbody></table>

The mass flow rate is calculated from the volumetric flow rate $\dot{V}$:

$$\dot{m} = \rho \dot{V} \qquad (48)$$

where $\rho$ is the fluid density. The volumetric flow rate $\dot{V}$ is obtained from hydraulic or aeraulic calculations. The infiltration (or the airtightness) of the building is measured by blower door. It may expressed in "air changes per hour" which represents the number of volumes of air contained by the building which are vehiculated in an hour. Some indicative values are given in Table 3. Typical values for airflow per person are given in

Table 2
<table><caption>Table 2 Typical values of airflow rate per person (Recknagel, et al., 2007) §3.5.1.1.2</caption><thead><tr><th>Typical situations</th><th>Volumetric airflow<br>m³/h per person</th></tr></thead><tbody><tr><td>Theatre, concert, cinema, exhibition halls, supermarkets, museums,<br>gyms</td><td>20</td></tr><tr><td>Restaurants, rest area, conference rooms, classrooms, auditorium</td><td>30</td></tr><tr><td>Office rooms</td><td>40</td></tr><tr><td>Open office rooms</td><td>60</td></tr></tbody></table>

## 2.3 Heat capacities

The heat capacity $C$ of a control volume is the amount of heat $Q$ that needs to be added in order to increase the temperature of $\Delta\theta$:

$$C = \lim_{\Delta\theta \to 0} \frac{\Delta Q}{\Delta \theta} \qquad (49)$$

For a homogeneous object of mass $m$ having the specific heat $c$,

$$C = mc \qquad (50)$$

From equation (49), the flow rate entering the body is

$$q = C\dot{\theta} \qquad (51)$$

## 3 Assembling of thermal networks

Buildings are systems composed of elements such as walls, windows, doors, etc. connected through heat and mass transfer. Therefore, the models of whole buildings are obtained by combining the models of individual components. Two important methods for obtaining large models are coupling and assembling. In coupling, the system of equations, which is obtained from the models of each element, is solved iteratively (Clarke, et al., 2007; El Khoury, et al., 2013; Ma
& Le, 1992). For example, building energy software tools use iterative methods to solve the coupled equations: EnergyPlus uses Gauss-Seidel successive substitution or Newton-Raphson method (Ficher, et al., 1999; EnergyPlus, 2015), TRNSYS uses successive substitutions and Powell's method (TRNSYS 17, 2009), ESP-r solves independently the domain equations and then the coupling (Clarke, et al., 2007), IDA ICE uses a modular approach (Sahlin, 1996).

The assembling of models is very different of coupling. While in coupling the models of the elements are separate, in assembling the complex model is a system of linear equations, at least at a certain moment during an iterative solving procedure. Assembling is an important reason for the use of models such as transfer functions in thermal modelling of buildings. More generally, input-output linear time invariant (LTI) models, such as state-space, transfer function, zero-pole-gain models (MathWorks, 2017) or two-port networks (Piotrowska & Chochowski, 2013; Ghosh, 2005; Chen, et al., 2013; Maillet, et al., 2000) may be connected to obtain a new, more complex, model. The model obtained by assembling has the advantage that can be analysed (e.g. find the eigenvalue and the time constants, the static gain, stability, controllability, observability, identifiability, etc.). However, these techniques are not applicable to networks or circuits that model transport phenomena in which the connections are done by conservation laws (such as conservation of mass, energy, momentum, electrical charge, etc.). The usual technique used for network models is coupling.

Circuits, networks, or bond graphs are widely used for modelling transfer phenomena (Karnopp, et al., 1990; Strang, 2007). The method of thermal networks (or circuits) is present in almost any primer on heat transfer. The heat conduction equation, introduced by Fourier, has been used for about two centuries to describe diffusion phenomena in dynamical physical systems. Ohm work on electricity was inspired by Fourier's heat conduction model; he considered the flow of electricity as being exactly analogue to the flow of heat. Fick also used an analogy with Fourier equation for transient diffusion of solutes in liquids. Models influenced by the diffusion equation are used for diffusion of gases, Brownian motion, flow in porous materials, random walk, etc. (Narasimhan, 1999). Therefore, it is important to have a procedure for assembling the networks (or the circuits).

There are well established algorithms for assembling models represented by finite elements (Nikishkov, 2010; Cueva-Zepeda & Avalos-Garcia, 2006; Ramabathiran & Gopalakrishnan, 2014).
However, for thermal networks the solution is obtained by solving iteratively the set of equations of the elementary models (Ruiz-Calvo, et al., 2016).

This section proposes a data structure for thermal networks and an assembling procedure. Then, aspects of software implementation are discussed. The example is given for a very simple, yet relevant, network on which the procedure can be checked by hand.

## 3.1 Defining the problem of circuit assembling

Given a number of thermal circuits, $TC_1, TC_2, ..., TC_n$ and knowing that some of their nodes are in common, find the assembled circuit $TC$. A simple example is given in Figure 5. There are four thermal circuits, $TC_1, TC_2, TC_3, TC_4$, having in common some of their nodes: the node 5 of $TC_1$ is common with the node 1 of $TC_2$, the node 2 of $TC_2$ is common with node 2 of $TC_3$ and the node 3 of $TC_2$ is in common with the node 2 of $TC_4$ (Figure 5b). Find the model of the assembled circuit $TC$ shown in Figure 5a.

From conservation of energy, it results that if there is a flow source in the node of the assembled circuit $TC$, it needs to be the sum of the sources in the respective nodes of each circuit $TC_k$. For example, the flow source in the node 5 of the assembled circuit from Figure 5a is the sum of flow sources present in node 5 of $TC_1$ and the flow source present in node 1 of $TC_2$. Since the thermal capacity is proportional to mass, from the conservation of mass, it results that if there is a capacity in a node of the assembled circuit, it needs to be the sum of
<figure><img src="image_9.png" /></figure>

<figure><img src="image_10.png" /><figcaption>Assembly-connectivity matrix</figcaption></figure>

$$ \mathbf{A}_c = \begin{bmatrix} 1 & 5 & 2 & 1 \\ 2 & 2 & 3 & 2 \\ 2 & 3 & 4 & 1 \end{bmatrix} $$

<figure><img src="image_11.png" /><figcaption>Figure 5 Example of the problem of assembling thermal circuits: given four circuits, assemble them knowing the common nodes. a) Assembled circuit. b) Four disassembled circuits. c) Assembling matrix. d) Algebraic description of each disassembled circuit.</figcaption></figure>
the capacities in the respective nodes of each circuit TC<sub>k</sub>. For example, the capacity in node 5 of the assembled circuit from Figure 5a is the sum of capacities present in node 5 of TC<sub>1</sub> and in node 1 of TC<sub>2</sub> shown in Figure 5b.

To exemplify the procedure, we will use a toy model representing a building formed by an insulated concrete wall and a glass wall. The room is ventilated and its air temperature is controlled by a P-controller. Auxiliary load is added to the room (Figure 6). The toy model is used to show specific aspects of the assembling procedure, not for the correctness of the modelling.

<figure><img src="image_12.png" /><figcaption>Figure 6 Toy model used for example</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 7 Model for components to be assembled</figcaption></figure>
We would like to construct separate models for concrete wall, glass wall, ventilation, and room air (Figure 7) and to assemble them into one model (Figure 6).

## 3.2 Algebraic description of the thermal circuits

A circuit is a weighted oriented graph with node representing temperatures, branches representing heat flows, and the weights representing the thermal conductances. Some nodes have thermal capacities, but not all of them. Some branches have temperature sources and some nodes have flow sources, but not all of them. The sources represent the inputs (i.e. the independent) variables of the model. The temperatures of the nodes and the flows in the branches represent the unknowns for which the problem is solved. Usually, only some of the temperatures of the nodes and/or flow in the branches are of practical interest and represent the output of the model. If the heat flow rate of a branch is considered as an output, then the temperatures in the nodes of the branch need to be found; the flow in the branch is calculated as the product between the conductance and the difference of temperatures.

A thermal circuit may be described by three matrices and three vectors. The matrices are:

1) **A** is an oriented incidence matrix with the number of rows equal to the number of branches and the number of columns equal to the number of nodes of the thermal circuit. The elements of matrix **A** are:

$$ a_{ij} = \begin{cases} 0 & \text{if the heat flow rate } i \text{ is not connected to the node } j \\ -1 & \text{if the heat flow rate } i \text{ leaves the node } j \\ 1 & \text{if the heat flow rate } i \text{ enters the node } j \end{cases} \qquad (52) $$

2) **G** is a diagonal matrix of conductances of dimension equal to the number of rows of **A**, i.e. the number of branches or the number of conductances. The elements of matrix **G** are:

$$ g_{ij} = \begin{cases} R_i^{-1} & \text{for } i = j \\ 0 & \text{for } i \neq j \end{cases} \qquad (53) $$

Note that each branch needs to have a conductance.

3) **C** is a diagonal matrix of capacitances of dimension equal to the number of columns of **A**, i.e. the number of nodes of the thermal circuit. The elements of matrix **C** are:

$$ c_{ij} = \begin{cases} C_i & \text{for } i = j \\ 0 & \text{for } i \neq j \end{cases} \qquad (54) $$
Note that not all nodes have a thermal capacity. Since the thermal capacity is proportional to mass and the mass proportional to volume, a node representing a surface will always have a zero capacity.

The vectors are:

1) **b** is a vector indicating the branches which have temperature sources. Its size is equal to the number of rows of matrix **A**, i.e. the number of branches. Its elements are:

$$b_i = \begin{cases} 1 & \text{if there is a temperature source on branch } i \\ 0 & \text{otherwise} \end{cases} \qquad (55)$$

2) **f** is a vector indicating the nodes which have a heat flow rate sources. Its size is equal to the number of columns of matrix **A**, i.e. the number of nodes. Its elements are:

$$f_i = \begin{cases} 1 & \text{for flow source in node } i \\ 0 & \text{otherwise} \end{cases} \qquad (56)$$

3) **y** is a vector indicating the temperature that are considered as outputs. Its size is equal to vector **f**. Its elements are:

$$y_i = \begin{cases} 1 & \text{for temperature of node } i \text{ as output variable} \\ 0 & \text{otherwise} \end{cases} \qquad (57)$$

Any thermal circuit **TC** can be described by the list of arrays: {**A**, **G**, **C**, **b**, **f**, **y**} (Figure 5d).

## 3.3 Numbering the thermal circuits

The construction of the list of arrays {**A**, **G**, **C**, **b**, **f**, **y**} requires the numbering of circuits. In principle, the numbering of the nodes and branches can be done arbitrarily. Once the numbering of the elementary circuits is done, the numbering of the assembled circuit is automatic.

### 3.3.1 Numbering elementary circuits

The connections between nodes are indicated by the oriented incidence matrix **A**. Since numbering becomes tedious for large circuits, the following rules may be adopted (Figure 5b):

- number the nodes in order (from left to right or from right to left;

- number the branches in increasing order of nodes and orient them from the lower to the higher node. Note: reference temperature is node 0.
As an example, for the thermal circuit $TC_1$ (red in Figure 5b), the branches are:

<table><thead><tr><th>1<sup>st</sup> node</th><th>2<sup>nd</sup> node</th><th>Branch</th></tr></thead><tbody><tr><td>0</td><td>1</td><td>1</td></tr><tr><td>1</td><td>2</td><td>2</td></tr><tr><td>2</td><td>3</td><td>3</td></tr><tr><td>3</td><td>4</td><td>4</td></tr><tr><td>4</td><td>5</td><td>5</td></tr></tbody></table>

For $TC_2$ (blue in Figure 5b), the branches are:

<table><thead><tr><th>1<sup>st</sup> node</th><th>2<sup>nd</sup> node</th><th>Branch</th></tr></thead><tbody><tr><td>1</td><td>2</td><td>1</td></tr><tr><td>1</td><td>3</td><td>2</td></tr><tr><td>2</td><td>3</td><td>3</td></tr></tbody></table>

For $TC_3$ (green in Figure 5b the branches are:

<table><thead><tr><th>1<sup>st</sup> node</th><th>2<sup>nd</sup> node</th><th>Branch</th></tr></thead><tbody><tr><td>0</td><td>1</td><td>1</td></tr><tr><td>1</td><td>2</td><td>2</td></tr></tbody></table>

For $TC_4$ (violet in Figure 5b), the nodes have the same numbers:

<table><caption>Table 3 Local and global indexing of nodes</caption><thead><tr><th>1<sup>st</sup> node</th><th>2<sup>nd</sup> node</th><th>Branch</th></tr></thead><tbody><tr><td>0</td><td>1</td><td>1</td></tr><tr><td>0</td><td>1</td><td>2</td></tr></tbody></table>

<table><thead><tr><th>Thermal circuit</th><th>TC1</th><th>TC2</th><th>TC3</th><th>TC4</th></tr></thead><tbody><tr><td>Local node index</td><td>1 2 3 4 5</td><td>1 2 3</td><td>1 2</td><td>1</td></tr><tr><td>Global node index</td><td>1 2 3 4 5</td><td>5 6 7</td><td>8 6</td><td>7</td></tr></tbody></table>
### 3.3.2 Numbering the assembled circuit

When assembling the thermal circuits, some nodes are put in common. Therefore, the number of nodes in the assembled circuit will be smaller than the sum of the nodes of elementary circuits. The number of branches will not change. The nodes and the branches of the assembled circuit will be in the order of assembling (Figure 5a, Table 3).

The assembling of the circuits is indicated by the assembling matrix. Each row of this matrix has four elements that indicate two nodes that will be put together:

- number of circuit 1,

- node of circuit 1,

- number of circuit 2,

- node of circuit 2.

For our example, the assembling matrix is:

$$ \mathbf{Ass} = \begin{bmatrix} 1 & 5 & 2 & 1 \\ 2 & 2 & 3 & 2 \\ 2 & 3 & 4 & 1 \end{bmatrix} \qquad (58) $$

The description of the assembled circuit, given by the list $TC = \{TC_1, ..., TC_i\}$ of list of arrays $TC_i = \{A_i, G_i, C_i, b_i, f_i, y_i\}$ (Figure 5c), and the assembling matrix **Ass** contain all the necessary information for obtaining the assembled circuit.

## 3.4 Assembling the circuits

The analysis (or the direct problem) of a thermal circuit $TC_i$ is to solve for $q_i$ and $\theta_i$ the equation:

$$ \begin{bmatrix} \mathbf{G}_i^{-1} & \mathbf{A}_i \\ -\mathbf{A}_i^T & \mathbf{C}_i s \end{bmatrix} \begin{bmatrix} \mathbf{q}_i \\ \mathbf{\theta}_i \end{bmatrix} = \begin{bmatrix} \mathbf{b}_i \\ \mathbf{f}_i \end{bmatrix} \qquad (59) $$

or find $\mathbf{u}_i$ from equation:

$$ \mathbf{K}_i \mathbf{u}_i = \mathbf{a}_i \qquad (60) $$

where

$$ \mathbf{K}_i = \begin{bmatrix} \mathbf{G}_i^{-1} & \mathbf{A}_i \\ -\mathbf{A}_i^T & \mathbf{C}_i s \end{bmatrix}; \mathbf{u}_i = \begin{bmatrix} \mathbf{q}_i \\ \mathbf{\theta}_i \end{bmatrix}; \mathbf{a}_i = \begin{bmatrix} \mathbf{b}_i \\ \mathbf{f}_i \end{bmatrix}. $$

Let's note the disassembled block vectors $\mathbf{u}_d, \mathbf{a}_d$ and matrix $\mathbf{K}_d$:
$$ \mathbf{u}_d = \begin{bmatrix} \mathbf{u}_1 \\ \vdots \\ \mathbf{u}_n \end{bmatrix}; \quad \mathbf{a}_d = \begin{bmatrix} \mathbf{a}_1 \\ \vdots \\ \mathbf{a}_n \end{bmatrix}; \quad \mathbf{K}_d = \begin{bmatrix} \mathbf{K}_1 & 0 & 0 \\ 0 & \ddots & 0 \\ 0 & 0 & \mathbf{K}_n \end{bmatrix} $$

There is a disassembling matrix $\mathbf{A}_d$ which transforms the assembled vectors (i.e. the block vector of elementary circuits) into disassembled vectors:

$$ \mathbf{u}_d = \mathbf{A}_d \mathbf{u}; \quad \mathbf{a}_d = \mathbf{A}_d \mathbf{a}; \qquad (61) $$

The assembled matrix and vectors are obtained by using disassembling matrix $\mathbf{A}_d$:

$$ \mathbf{K} = \mathbf{A}_d^T \mathbf{K}_d \mathbf{A}_d \qquad (62) $$

$$ \mathbf{u} = \mathbf{A}_d^T \mathbf{u}_d \qquad (63) $$

$$ \mathbf{a} = \mathbf{A}_d^T \mathbf{a}_d \qquad (64) $$

The elements of the assembled circuit, $\mathbf{A}, \mathbf{G}, \mathbf{C}, \mathbf{b}, \mathbf{f}, \mathbf{y}$, are then obtained from:

$$ \mathbf{K} = \begin{bmatrix} \mathbf{G}^{-1} & \mathbf{A} \\ -\mathbf{A}^T & \mathbf{C}_S \end{bmatrix}; \mathbf{u} = \begin{bmatrix} \mathbf{q} \\ \mathbf{\theta} \end{bmatrix}; \mathbf{a} = \begin{bmatrix} \mathbf{b} \\ \mathbf{f} \end{bmatrix} \qquad (65) $$

## 3.5 Algorithm

### 3.5.1 Obtaining the global indexes of the assembling matrix

In order to indicate the common nodes of the circuits, it is convenient to give the assembling matrix $\mathbf{Ass}$, of which an example is given in equation (58), with 4 elements on each line:

1. number of the 1st circuit,

2. local number of the node of the 1st circuit,

3. number of the 2nd circuit,

4. local number of the node of the 2nd circuit.
<figure><img src="image_14.png" /><figcaption>Figure 8 The disassembling matrix: rows correspond to disassembled circuits; columns correspond to assembled circuits.</figcaption></figure>

We need to obtain an assembling matrix **A**<sub>ss</sub> of two columns of global disassembled nodes that are put in common. For our example (Figure 5):

- the node 5 of $TC_1$ is put in common with the node 1 of $TC_2$, which has the global value 5 + 1 = 6 (5 = number of nodes of $TC_1$, 1 = local index in $TC_2$);

- the node 2 of $TC_2$ (global value 5 + 2) is put in common with the node 2 of $TC_3$ (global value 5 + 3 + 2 = 10, where 5 = number of nodes of $TC_1$, 3 = number of the nodes of $TC_2$, 2 local index in $TC_3$);
- the node 3 of $TC_2$ (global value $5+3=8$) is put in common with the node 1 of $TC_4$ (global value $5+3+2+1=11$, where 5 = number of nodes of $TC_1$, 3 = number of the nodes of $TC_2$, 2 = number of nodes in $TC_3$, 1 = local index in $TC_4$);

From

$$ \mathbf{AssX} = \begin{bmatrix} 1 & 5 & 2 & 1 \\ 2 & 2 & 3 & 2 \\ 2 & 3 & 4 & 1 \end{bmatrix} \qquad (66) $$

we obtain:

$$ \mathbf{Ass} = \begin{bmatrix} 5 & 6 \\ 7 & 10 \\ 8 & 11 \end{bmatrix} \qquad (67) $$

The information on the number of branches and nodes for each thermal circuit $TC_k$ is taken from the length of vector **b**.

### 3.5.2 Obtaining the disassembling matrix

There disassembling matrix $A_d$ transforms the assembled vectors into dissembled vectors (i.e. the block vector of elementary circuits):

$$ \mathbf{u}_d = A_d \mathbf{u}; \quad \mathbf{a}_d = A_d \mathbf{a}; \qquad (68) $$

The assembling implies that some of the nodes are merged: their number decreases and their "global" index changes.

First, create a block matrix that keeps the indexes of the temperature nodes:

$$ A_{d\theta} = \begin{bmatrix} \begin{bmatrix} 1 & & & \\ & 1 & & \\ & & 1 & \\ & & & 1 \end{bmatrix} & & & \\ & & & \\ & & & \\ & & & \end{bmatrix} \begin{bmatrix} 1 & & & \\ & 1 & & \\ & & 1 & \\ & & & 1 \end{bmatrix} \begin{bmatrix} 1 & & & \\ & 1 & & \\ & & 1 & \\ & & & 1 \end{bmatrix} \begin{bmatrix} 1 & & & \\ & 1 & & \\ & & 1 & \\ & & & 1 \end{bmatrix} \begin{bmatrix} 1 & & & \\ & 1 & & \\ & & 1 & \\ & & & 1 \end{bmatrix} \end{bmatrix} $$

Then add the columns that merge:
<figure><img src="image_15.png" /></figure>

Eliminate the columns that correspond to the eliminated nodes to obtain

<figure><img src="image_16.png" /></figure>

The branches (flows) keep their global number. The disassembling matrix is then obtained from the block matrix (Figure 9 a) by re-arranging the rows in order correspond to the vector (Figure 9 b): $[q_1^T \quad \theta_1^T \quad q_2^T \quad \theta_2^T \quad q_3^T \quad \theta_3^T \quad q_4^T \quad \theta_4^T]^T$

### 3.5.3 Algorithm for the disassembling matrix

Having the disassembling matrix $A_d$, the assembling

$$ K = A_d^T K_d A_d \tag{69} $$

needs the matrix $K_d$ which is a block matrix

$$ K_d = \begin{bmatrix} K_1 & 0 & 0 \\ 0 & \ddots & 0 \\ 0 & 0 & K_n \end{bmatrix} $$

of block matrices of each thermal circuit $TC_i$

$$ K_i = \begin{bmatrix} G_i^{-1} & A_i \\ -A_i^T & C_i s \end{bmatrix} $$
<figure><img src="image_17.png" /><figcaption>Figure 9 Obtaining the disassembling matrix:</figcaption></figure>

a) block matrix for $[\mathbf{q}_1^T \quad \mathbf{q}_2^T \quad \mathbf{q}_3^T \quad \mathbf{q}_4^T \quad \boldsymbol{\theta}_1^T \quad \boldsymbol{\theta}_2^T \quad \boldsymbol{\theta}_3^T \quad \boldsymbol{\theta}_4^T]^T$

b) block matrix re-arranged for $[\mathbf{q}_1^T \quad \boldsymbol{\theta}_1^T \quad \mathbf{q}_2^T \quad \boldsymbol{\theta}_2^T \quad \mathbf{q}_3^T \quad \boldsymbol{\theta}_3^T \quad \mathbf{q}_4^T \quad \boldsymbol{\theta}_4^T]^T$
<figure><img src="image_18.png" /><figcaption>Algorithm for obtaining the dissembling matrix</figcaption></figure>

Figure 10 Obtain matrix temperature nodes $A_{d\theta}$: a) create diagonal matrix; b) add the column of the 2^nd^ node to the column of the 1^st^ node in merging: [1 5 2 1] $\rightarrow$ col5=col5+col6; [2 2 3 2] $\rightarrow$ col7 = col7 + col10; [2 3 4 1] $\rightarrow$ col8 = col8 = col11; c) delete the columns of the 2^nd^ node in merging; 4) obtain matrix $A_{d\theta}$

<figure><img src="image_19.png" /><figcaption>Figure 11 Obtain matrix for flow branches $A_{dq}$</figcaption></figure>
<figure><img src="image_20.png" /><figcaption>Figure 12 Block matrix for [$q_1^T \ q_2^T \ q_3^T \ q_4^T \ \theta_1^T \ \theta_2^T \ \theta_3^T \ \theta_4^T$]$^T$ (first flow branches, then temperature nodes)</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 13 Matrix $A_d$ for [$q_1^T \ \theta_1^T \ q_2^T \ \theta_2^T \ q_3^T \ \theta_3^T \ q_4^T \ \theta_4^T$]$^T$ after rearranging in order to obtain in the order of thermal circuits</figcaption></figure>
### 3.5.4 Assembled circuits

The assembled matrix and vector are obtained by using the disassembling matrix $A_d$:

$$K = A_d^T K_d A_d \qquad (70)$$

The elements of the assembled circuit, A, G, C, b, f, y, are then obtained from:

$$K = \begin{bmatrix} G^{-1} & A \\ -A^T & C_S \end{bmatrix}; u = \begin{bmatrix} q \\ \theta \end{bmatrix}; a = \begin{bmatrix} b \\ f \end{bmatrix} \qquad (71)$$

### 3.5.5 Global assembled indexes

The global indexes of the assembled circuit (Figure 13) result from the calculation of the disassembling matrix $A_d$ (Figure 14).

<figure><img src="image_22.png" /><figcaption>Figure 14 The relation between the local indexes of the blocks of circuits and the global indexes of the assembled circuit</figcaption></figure>

## 4 Transforming thermal circuits into state space representation

Thermal circuits are linear models with constant coefficients. However, in control theory, the state-space representation is widely used. The aim of this chapter is to transform a thermal circuit, formed by resistors, capacities, temperature sources and heat flow sources, into its state space representation. The problem in this transformation is that some capacities may be zero.

In this section we will use a very simple model of heat transfer through a wall (Figure 15). The model is very simplified (e.g. the number of meshes in the wall are too small) in order to keep the presentation manageable by hand calculations.
<figure><img src="image_23.png" /></figure>

**Legend**

*Temperature source*

$T_o$ outdoor temperature

*Heat flow sources*

*Q*<sub>o</sub> absorbed incident solar and long wave radiation on the outdoor surface;

*Q*<sub>i</sub> short wave radiation from lights, transmitted solar, long wave radiation exchange with other surfaces, and long wave radiation from internal sources on the indoor surface;

*Q*<sub>g</sub> heat flow gained by convection from the internal sources;

*Q*<sub>HVAC</sub> heat flow from the HVAC system, i.e. the thermal load.

*Temperatures in the thermal network*

$\theta_{os}, \theta_{is}$ outdoor and indoor surface temperatures, respectively,

$\theta_{w1}, \theta_{w2}$ temperatures in the wall.

*Heat flows in the thermal network*

$q_v$ heat flow rate through the window (advection, conduction, convection);

$q_{co}, q_{ci}$ convection from outside and inside air, respectively;

$q_{w1}, q_{w2}, q_{w3}$ conduction through the wall.

<caption>Figure 15 Typical thermal circuit for heat balance method: a) usual representation; b) representation in the form of typical branches</caption>
<figure><img src="image_24.png" /><figcaption>Figure 16 Temperature nodes are ordered so that the capacity matrix contains zero-blocks</figcaption></figure>

<figure><img src="image_25.png" /><figcaption>Figure 17 Temperature nodes need to be re-ordered</figcaption></figure>

<table><caption>Table 4 Parameter of the thermal network shown in Figure 1</caption><thead><tr><td>Parameter</td><td>Value</td></tr></thead><tbody><tr><td>Indoor air capacity, C<sub>a</sub></td><td>82 · 10<sup>3</sup> J/K</td></tr><tr><td>Wall capacity, C<sub>w1</sub> = C<sub>w2</sub></td><td>2 · 10<sup>6</sup> J/K</td></tr><tr><td>Thermal conductance of the wall, R<sub>w</sub><sup>-1</sup></td><td>1.45 W/K</td></tr><tr><td>Thermal conductance of one third of the wall, R<sub>w1</sub><sup>-1</sup> = R<sub>w2</sub><sup>-1</sup> = R<sub>w3</sub><sup>-1</sup> = 3 R<sub>w</sub><sup>-1</sup></td><td>4.35 W/K</td></tr><tr><td>Thermal conductance of the window and due to losses by ventilation, R<sub>v</sub><sup>-1</sup></td><td>38.3 W/K</td></tr><tr><td>Outdoor convection conductance, R<sub>co</sub><sup>-1</sup></td><td>250.0 W/K</td></tr><tr><td>Indoor convection conductance, R<sub>ci</sub><sup>-1</sup></td><td>125.0 W/K</td></tr></tbody></table>
$\theta_{so}$ $\theta_{si}$ $\theta_{w1}$ $\theta_{w2}$ $\theta_{a}$

$$
\begin{array}{c}
\begin{array}{c}
\begin{array}{c}
\begin{array}{c}
\theta_{so} & \theta_{si} & \theta_{w1} & \theta_{w2} & \theta_{a}
\end{array}
\\[1em]
\mathbf{A} = \begin{bmatrix}
1 & 0 & 0 & 0 & 0 \\
-1 & 0 & 1 & 0 & 0 \\
0 & 0 & -1 & 1 & 0 \\
0 & 1 & 0 & -1 & 0 \\
0 & -1 & 0 & 0 & 1 \\
0 & 0 & 0 & 0 & 1
\end{bmatrix}; \mathbf{G} = \begin{bmatrix}
R_{co}^{-1} & 0 & 0 & 0 & 0 \\
0 & R_{w1}^{-1} & 0 & 0 & 0 \\
0 & 0 & R_{w2}^{-1} & 0 & 0 \\
0 & 0 & 0 & R_{w3}^{-1} & 0 \\
0 & 0 & 0 & 0 & R_{ci}^{-1} \\
0 & 0 & 0 & 0 & R_{f}^{-1}
\end{bmatrix}; \mathbf{b} = \begin{bmatrix}
T_{ow} \\
0 \\
0 \\
0 \\
0 \\
T_{of}
\end{bmatrix}
\end{array}
\\[1em]
\begin{array}{c}
\mathbf{C} = \begin{bmatrix}
0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 \\
0 & 0 & C_{w1} & 0 & 0 \\
0 & 0 & 0 & C_{w2} & 0 \\
0 & 0 & 0 & 0 & C_{a}
\end{bmatrix} \\
\mathbf{f} = [ \begin{array}{c} \dot{Q}_o \\ 0 \\ 0 \\ \dot{Q}_i \\ \dot{Q}_a \end{array} ]^T
\end{array}
\\[1em]
\begin{array}{l}
\left\{
\begin{array}{l}
\mathbf{e} = -\mathbf{A}\theta + \mathbf{b} \\
\mathbf{q} = \mathbf{G}\mathbf{e} \\
\mathbf{C}\dot{\theta} = \mathbf{A}^T\mathbf{q} + \mathbf{f}
\end{array}
\right. \\
\mathbf{C}\dot{\theta} = -\mathbf{A}^T\mathbf{G}\mathbf{A}\theta + \mathbf{A}^T\mathbf{G}\mathbf{b} + \mathbf{f}
\end{array}
\\[1em]
\dot{Q}_a = \dot{Q}_g + \dot{Q}_{HVAC}
\end{array}
$$

<caption>Figure 18. Obtaining the differential-algebraic equations for the circuit from Figure 17</caption>

$\theta_{so}$ $\theta_{si}$ $\theta_{w1}$ $\theta_{w2}$ $\theta_{a}$

$$
\begin{array}{c}
\begin{array}{c}
\begin{array}{c}
\theta_{so} & \theta_{si} & \theta_{w1} & \theta_{w2} & \theta_{a}
\end{array}
\\[1em]
\mathbf{A} = \begin{bmatrix}
1 & 0 & 0 & 0 & 0 \\
-1 & 0 & 1 & 0 & 0 \\
0 & 0 & -1 & 1 & 0 \\
0 & 1 & 0 & -1 & 0 \\
0 & -1 & 0 & 0 & 1 \\
0 & 0 & 0 & 0 & 1
\end{bmatrix}; \mathbf{G} = \begin{bmatrix}
R_{co}^{-1} & 0 & 0 & 0 & 0 \\
0 & R_{w1}^{-1} & 0 & 0 & 0 \\
0 & 0 & R_{w2}^{-1} & 0 & 0 \\
0 & 0 & 0 & R_{w3}^{-1} & 0 \\
0 & 0 & 0 & 0 & R_{ci}^{-1} \\
0 & 0 & 0 & 0 & R_{f}^{-1}
\end{bmatrix}; \mathbf{b} = \begin{bmatrix}
T_{ow} \\
0 \\
0 \\
0 \\
0 \\
T_{of}
\end{bmatrix}
\end{array}
\\[1em]
\begin{array}{c}
\mathbf{C} = \begin{bmatrix}
0 & 0 & 0 & 0 & 0 \\
0 & C_{w1} & 0 & 0 & 0 \\
0 & 0 & C_{w2} & 0 & 0 \\
0 & 0 & 0 & C_{w3} & 0 \\
0 & 0 & 0 & 0 & C_{a}
\end{bmatrix} \\
\mathbf{f} = [ \begin{array}{c} \dot{Q}_o \\ 0 \\ 0 \\ \dot{Q}_i \\ \dot{Q}_a \end{array} ]^T
\end{array}
\\[1em]
\begin{array}{l}
\left\{
\begin{array}{l}
\mathbf{e} = -\mathbf{A}\theta + \mathbf{b} \\
\mathbf{q} = \mathbf{G}\mathbf{e} \\
\mathbf{C}\dot{\theta} = \mathbf{A}^T\mathbf{q} + \mathbf{f}
\end{array}
\right. \\
\mathbf{C}\dot{\theta} = -\mathbf{A}^T\mathbf{G}\mathbf{A}\theta + \mathbf{A}^T\mathbf{G}\mathbf{b} + \mathbf{f}
\end{array}
\\[1em]
\dot{Q}_a = \dot{Q}_g + \dot{Q}_{HVAC}
\end{array}
$$

<caption>Figure 19. Obtaining the differential-algebraic equations for the circuit from Figure 17</caption>

## 4.1 Obtaining the system of differential-algebraic equations

The system of equations corresponding to a thermal circuit may be obtained by using the Kirchhoff's laws and the constitutive laws for thermal transfer (Ghiaus, 2013; Strang,
2007). The steps needed to obtain the differential-algebraic system of equations may be
synthesized in the following algorithm:

1. Group the temperatures according to the type of node: without capacitance ($\theta_{so}$, $\theta_{si}$, $\theta_a$), and with capacities ($\theta_{w1}$, $\theta_{w2}$):

$$ \mathbf{\Theta} = [\theta_{so} \ \theta_{si} \ \theta_a \ \theta_{w1} \ \theta_{w2}]^T \qquad (72) $$

2. Write the matrices describing the circuit (see an example in Figure 18 and Figure 19 for the thermal circuit given in Figure 15):

a. Transform the thermal circuit in an oriented graph by indicating the direction of the heat transfer rate for each branch. The directions are arbitrary (if the nodes are numbered, can be in increasing order of the node numbering).

b. Write the oriented incidence matrix **A**, equation (52), conductance matrix **G**, equation (53), and capacitance matrix **C**, equation (54).

c. Write the vectors of temperature sources **b** given by equation (55), of flow rate sources **f** given by equation (56) and of outputs **y** given by equation (57).

3. Apply Kirchhoff's laws and the constitutive laws to obtain the differential-algebraic system of equations:

$$ C \dot{\mathbf{\Theta}} = -A^T G A \mathbf{\Theta} + A^T G b + f $$

By writing the differences of temperature according to equation (28) (equivalent to the Kirchhoff's voltage law),

$$ e = -A \mathbf{\Theta} + b \qquad (73) $$

the balance of heat rates in nodes (equivalent to the Kirchhoff's current law),

$$ C \dot{\mathbf{\Theta}} = A^T \mathbf{q} + f \qquad (74) $$

and the constitutive laws for heat transfer,

$$ \mathbf{q} = G e \qquad (75) $$

we obtain the differential algebraic equations describing the thermal circuit:

$$ C \dot{\mathbf{\Theta}} = K \mathbf{\Theta} + K_b \mathbf{b} + f; \qquad (76) $$

where

$$ K = -A^T G A \mathbf{\Theta} \text{ and } K_b = A^T G \qquad (77) $$
If the diagonal matrix C has elements on the diagonal which are zero, the system of equations (76) is a system of differential algebraic equations.

## 4.2 Obtaining the state space representation from the thermal circuit

If the thermal circuit contains nodes without capacitance, the matrix C is singular. In order to obtain the state space model, the equations corresponding to the nodes without capacitance need to be eliminated from the system of equations (76) (Ghiaus, 2013). By partitioning the matrix C,

$$C = \begin{bmatrix} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & c_C \end{bmatrix} \qquad (78)$$

where $c_C$ corresponds to the nodes having capacities, the set of equations (76) may be written as:

$$\begin{bmatrix} \mathbf{0} & \mathbf{0} \\ \mathbf{0} & c_C \end{bmatrix} \begin{bmatrix} \dot{\theta}_0 \\ \dot{\theta}_C \end{bmatrix} = \begin{bmatrix} \mathbf{K}_{11} & \mathbf{K}_{12} \\ \mathbf{K}_{21} & \mathbf{K}_{22} \end{bmatrix} \begin{bmatrix} \theta_0 \\ \theta_C \end{bmatrix} + \begin{bmatrix} \mathbf{K}_{b1} \\ \mathbf{K}_{b2} \end{bmatrix} \mathbf{b} + \begin{bmatrix} \mathbf{I}_{11} & \mathbf{0} \\ \mathbf{0} & \mathbf{I}_{22} \end{bmatrix} \begin{bmatrix} \mathbf{f}_0 \\ \mathbf{f}_C \end{bmatrix} \qquad (79)$$

where

$\theta_0$ and $f_0$ correspond to the nodes without thermal capacity;

$\theta_C$ and $f_C$ correspond to the nodes with thermal capacity;

$C_C$ is the bloc of the partitioned matrix C for which the elements on the diagonal are non-zero;

$K_{11}$, $K_{12}$, $K_{21}$, and $K_{22}$ are blocs of the partitioned matrix K obtained according to the partitioning of the matrix C;

$K_{b1}$ and $K_{b2}$ are blocs of the partitioned matrix $K_b$ obtained according to the partitioning of the matrix C;

$I_{11}$ and $I_{22}$ are identity matrices.

The state equation of the state space model is

$$\dot{\theta}_C = A_S \theta_C + B_S u \qquad (80)$$

where the state matrix is

$$A_S = C_C^{-1}(-K_{21}K_{11}^{-1}K_{12} + K_{22}) \qquad (81)$$

and the input matrix is
$$ B_s = C_c^{-1} [-K_{21}K_{11}^{-1}K_{b1} + K_{b2} - K_{21}K_{11}^{-1}I] \quad (82) $$

For the numerical example of the thermal circuit shown in Figure 15 with the values of its parameters given in Table 4, the state variables are

$$ \theta = [\theta_a \ \theta_{w1} \ \theta_{w2}]^T \quad (83) $$

and the bloc vector of inputs is:

$$ u = [b \ f_0 \ f_c]^T \quad (84) $$

The numerical values of the matrices of the model are:

$$ C = \begin{bmatrix} 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 \\ 0 & 82 \cdot 10^3 & 0 & 0 & 0 \\ 0 & 0 & 2 \cdot 10^6 & 0 & 0 \\ 0 & 0 & 0 & 2 \cdot 10^6 \end{bmatrix} = \begin{bmatrix} 0 & 0 \\ 0 & C_c \end{bmatrix} \quad (85) $$

$$ K = \begin{bmatrix} -254.35 & 0 & 0 & 4.35 & 0 \\ 0 & -129.35 & 125 & 0 & 4.35 \\ 0 & 125 & -163.30 & 0 & 0 \\ 4.35 & 0 & 0 & -8.70 & 4.35 \\ 0 & 4.35 & 0 & 4.35 & -8.70 \end{bmatrix} = \begin{bmatrix} K_{11} & K_{12} \\ K_{21} & K_{22} \end{bmatrix} \quad (86) $$

$$ K_b = \begin{bmatrix} 250 & -4.35 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 4.35 & -125 & 0 \\ 0 & 0 & 0 & 0 & 125 & 38.3 \\ 0 & 4.35 & -4.35 & 0 & 0 & 0 \\ 0 & 0 & 4.35 & -4.35 & 0 & 0 \end{bmatrix} = \begin{bmatrix} K_{b1} \\ K_{b2} \end{bmatrix} \quad (87) $$

Substituting these matrices in equation (81), we obtain the state matrix:

$$ A_s = \begin{bmatrix} -5.18 \cdot 10^{-4} & 0 & 5.13 \cdot 10^{-5} \\ 0 & -4.31 \cdot 10^{-6} & 2.17 \cdot 10^{-5} \\ 2.10 \cdot 10^{-6} & 2.17 \cdot 10^{-5} & -4.28 \cdot 10^{-5} \end{bmatrix} \quad (88) $$

The input matrix $B_s$ has 11 columns, corresponding to the input vector

$$ u = [b^T \ f_0^T \ f_c^T]^T = [T_{ow} \ 0 \ 0 \ 0 \ 0 \ T_{ov} \ \dot{Q}_o \ \dot{Q}_i \ \dot{Q}_{aux} \ 0 \ 0] \quad (89) $$
The inputs corresponding to zeros are useless and can be eliminated. Therefore, from the eleven columns of the input matrix $\mathbf{B}_s$ it can be retained only the columns corresponding to inputs 1, 6, 7, 8, and 9:

$$ \mathbf{B}_s = \begin{bmatrix} 0 & 4.67 \cdot 10^{-4} & 0 & 1.18 \cdot 10^{-5} & 1.22 \cdot 10^{-5} \\ 2.14 \cdot 10^{-6} & 0 & 8.55 \cdot 10^{-9} & 0 & 0 \\ 0 & 0 & 0 & 1.68 \cdot 10^{-8} & 0 \end{bmatrix} \quad (90) $$

The output matrix of the state-space representation extracts the output, i.e. the indoor air temperature $\theta_a$ which is the first element from the state vector:

$$ \mathbf{C}_s = \begin{bmatrix} 1 & 0 & 0 \end{bmatrix} \quad (91) $$

For the thermal circuit given in Figure 15, the feed-through matrix is zero,

$$ \mathbf{D}_s = \mathbf{0} \quad (92) $$

The complete state-space representation of the thermal circuit from Figure 15, with the values of its parameters given in Table 4, is

$$ \begin{cases} \dot{\theta}_c = \mathbf{A}_s \theta_c + \mathbf{B}_s \theta_c \\ \theta_a = \mathbf{C}_s \theta_c + \mathbf{D}_s \mathbf{u} \end{cases} \quad (93) $$

with the values of the matrices given by equations (88), (90), and (91).

If the term $C_a = m_a c_a \dot{\theta}_a$ is zero, then the state vector is

$$ \dot{\mathbf{\theta}} = [\theta_{w1} \quad \theta_{w2}]^T \quad (94) $$

which implies that the matrices in equations (85) - (87) are partitioned correspondingly. In our numerical example, the first three equations need to be eliminated from the system of equations (76). The state matrix $\mathbf{A}_s$, obtained with the expression (81), is

$$ \mathbf{A}_s = \begin{bmatrix} -4.31 \cdot 10^{-6} & 2.17 \cdot 10^{-6} \\ 2.17 \cdot 10^{-6} & -4.07 \cdot 10^{-6} \end{bmatrix} \quad (95) $$

The input matrix $\mathbf{B}_s$, obtained with the expression (82) and retaining the only the columns corresponding to inputs 1, 6, 7, 8, and 9, is:

$$ \mathbf{B}_s = \begin{bmatrix} 2.14 \cdot 10^{-6} & 0 & 8.55 \cdot 10^{-9} & 0 & 0 \\ 0 & 1.89 \cdot 10^{-6} & 0 & 6.46 \cdot 10^{-8} & 4.95 \cdot 10^{-8} \end{bmatrix} \quad (96) $$
The observation equation can be obtained from the first row of equation (79) (Ghiaus, 2013)

$$ \theta_0 = -K_{11}^{-1}(K_{12}\theta_c + K_{b1}\mathbf{b} + I_{11}\mathbf{f}_0) = -K_{11}^{-1}\left(K_{12}\theta_c + [K_{b1} \quad I_{11} \quad 0]\begin{bmatrix} \mathbf{b} \\ \mathbf{f}_0 \\ \mathbf{f}_c \end{bmatrix}\right) \quad (97) $$

Then, the output equation is

$$ C_s = -K_{11}^{-1} K_{12} \quad (98) $$

and the feed through matrix is

$$ D_s = -K_{11}^{-1}[-K_{b1} \quad I_{11} \quad 0] \quad (99) $$

Keeping only the non-zero inputs [$T_{ow}$, $T_{ov}$, $\dot{Q}_o$, $\dot{Q}_i$, $\dot{Q}_{aux}$] from the input vector $\mathbf{u}$ given by equation (89), the output and feed through matrices for our numerical example are:

$$ C_s = [0 \quad 9.89 \cdot 10^{-2}] \quad (100) $$

and

$$ D_s = [0 \quad 9.01 \cdot 10^{-1} \quad 0 \quad 2.27 \cdot 10^{-2} \quad 2.35 \cdot 10^{-2}] \quad (101) $$

## 5 Conclusions

Thermal networks are widely used to model heat transfer. The phenomena of conduction, convection, radiation and advection can be linearized; as a consequence, the thermal networks, which are weighted oriented graphs, can be represented by matrices and vectors.

State space representation is widely used in the analysis and synthesis of control systems. Linear time invariant models may be used as local linearized models of non-linear
systems. Therefore, the linear algebra representation of state-space models has a large field of applications.

This chapter described succinctly the modelling of heat transfer by thermal networks and emphasized the fact that the temperature and heat-flow rate sources are inputs (or independent variables), the temperatures in nodes and the flow rates are outputs (or dependent variables), and the resistances and the capacities are parameters of a model structure. A data structure composed of matrices and vectors was proposed. The novelty of the data structure is the definition of an output vector which represents a subset of temperatures that are needed as observables.

A second novelty presented in this chapter is the assembling of thermal circuits. This technique allows us to construct large models from constitutive blocs. For example, the model of a complex building may be obtained by interconnecting typical blocs such as walls, floors, doors, windows, etc. Complex systems can be obtained also by coupling the equations of the typical blocs and solving iteratively the system of equations. The advantage of assembling is that the model of the whole system is a single thermal network that can be analysed. The key point in assembling is obtaining the disassembling matrix. An algorithm for obtaining it is presented.

The third novelty is the transformation of thermal circuit in state space representation. While examples for simple circuits are abundant and other methods (such as Kirchhoff laws, nodal analysis, etc.) are available, the method proposed is directly related to the matrix representation of the thermal circuits. The principal characteristic of the method is Gauss elimination of the block-matrices and vectors related to node temperatures that do not have capacities (i.e. that are not state variables).

Assembling thermal circuits and obtaining state-space models from them can be used in at least two important fields. The first is Building Information Modelling (BIM): each component has its model and the model of the building can be obtained by assembling the models of components. The second is system theory in which the state space is a suitable form of the model for system analysis and synthesis: eigen-value decomposition, model
order reduction, model predictive control, observability, controllability and identifiability, etc.

Acknowledgement

The author highly appreciates the discussions with Florent Alzetto, research engineer at Saint Gobain Research. His feed-back and comments were very valuable.

6 References

Anon., 2016. Von Neumann stability analysis. [Online]

Available at: https://en.wikipedia.org/wiki/Von_Neumann_stability_analysis

Anon., n.d. [Online].

ASHRAE, 2001 §30. Fenestration. In: ASHRAE Handbook: Fundamentals. Atlanta, GA, USA: s.n.

ASHRAE, 2001 §38. Physical properties of materials. In: ASHRAE Handbook: Fundamentals. Atlanta, GA, USA: s.n., p. 544.

Bergman, T., Lavine, A., Incropera, F. & Dewitt, D., 2011. Fundamentals of Heat and Mass Transfer. 7 ed. s.l.:John Wiley & Sons.

Chadha, H. & Baugh, J., 1996. Network-distributes finite element analysis. Advances in Engineering Software, Volume 25, pp. 267-280.

Chen, Y., Athienitis, A. & Galal, K., 2013. Frequency domain and finite difference modeling of ventilated concrete slabs and comparison with field measurements: Part 1, modeling methodology. International Journal of Heat and Mass Transfer, pp. 948-956.

Clarke, J. A., 2001. Energy simulation in Building Design. 2 ed. Oxford: Butterworth - Heinemann.

Clarke, J., Kelly, N. J. & Tang, D., 2007. A review of ESP-r flexible solution approaches and its application to prospective technical domain developments. Advances in Building Energy Research, 1(1), pp. 227-247.

Cueva-Zepeda, A. & Avalos-Garcia, J., 2006. The role of a connectivity matrix in the assemblage process of the finite element method. Advances in Engineering Software, 37(11), pp. 721-727.

Ebery, D., 2008. Stability Analysis for Systems of Differential Equations. [Online]

Available at: https://www.geometrictools.com/Documentation/StabilityAnalysis.pdf
El Khoury, K., Mouawad, G., El Hitti, G. & Nemer, M., 2013. The component interaction network approach for modeling of complex thermal systems. *Advances in Engineering Software*, Volume 65, pp. 149-157.

EnergyPlus, 2012. *Engineering Reference*. s.l.:University of Illinois and Lawrence Berkeley National Laboratory.

EnergyPlus, 2015. *Engineering Reference*. University of Illinois: University of Illinois & LBNL.

EnergyPlus, 2015. *Engineering Reference*. s.l.:DOE.

Ficher, D. E. et al., 1999. *A modular loop-based approach to HVAC energy simulation and its implementation in EnergyPlus*. Kyoto, Japan, Building Simulation, pp. B-31.

Ghiaus, C., 2013. Causality issue in the heat balance method for calculating the design heating and cooling load. *Energy*, Volume 50, pp. 292-301.

Ghosh, S., 2005. *Network Theory: Analysis and Synthesis*. New Delhi: Pretince Hall of India.

Karnopp, D. C., Margolis, D. L. & Rosenberg, R. C., 1990. *System dynamics: a unified approach*. s.l.:Wiley.

Kocak, S. & Akay, H., 2001. Parallel Schur complement method for large-scale systems on distributed memory computers. *Applied Mathematical Modelling*, 25(10), pp. 873-886.

Maillet, D. et al., 2000. *Thermal Quadrupoles: Solving the Heat Equation through Integral Transforms*. s.l.:Wiley.

Ma, J. & Le, M., 1992. A new method for coupling of boundary element method and finite element method. *Applied Mathematical Modelling*, 16(1), pp. 43-46.

Marcia, N., 2005. *Modeling and control of dynamic systems*. Clifton Park, NY: Thomson Delmar Learning.

MathWorks, 2017. *Model Interconnection*. [Online]

Available at: https://mathworks.com/help/control/model-interconnection.html

Narasimhan, T. N., 1999. Fourier's heat conduction equation: History, influence, and connections. *Reviews of Geophysics*, 37(doi.org/10.1029/1998RG900006), pp. 151-172.

Naveros, I., Ghiaus, C., Ordonez, J. & Ruiz, D. P., 2016. Thermal networks from heat equation by using the finite element method. *Heat Transfer*, Volume 106, pp. 33-43.

Nikishkov, G., 2010. *Programming Finite Elements in Java*. s.l.:Springer.

Piotrowska, E. & Chochowski, A., 2013. Representation of transient heat transfer as the equivalent thermal network (ETN). *International Journal of Heat and Mass Transfer*, Volume 63, pp. 113-119.
Ramabathiran, A. A. & Gopalakrishnan, S., 2014. Automatic finite element formulation and assembly of hyperelastic higher order structural models. *Applied Mathematical Modelling*, pp. 2867-2883.

Recknagel, H., Sprenger, E. & Schramek, E.-R., 2007. *Génie Clmatique*. 72 ed. Paris: Dnod.

Ruiz-Calvo, F. et al., 2016. Coupling short-term (B2G model) and long-term (g-function) models for ground source heat exchanger simulation in TRNSYS. Application in a real installation. *Applied Thermal Engineering*, Volume 102, pp. 720-732.

Sahlin, P., 1996. *Modelling and Simulation Methods for Modular Continuous Systems in Buildings*, Stockholm: Royal Institute of Technology.

Save, Y. D., Narayanan, H. & Patkar, S. B., 2011. Solution of Partial Differential Equations by electrical analogy. *Journal of Computational Science*, 2(1), pp. 18-30.

Strang, G., 2007. *Computational science and engineering*. Wellesley CA: Wellesley-Cambridge.

Strang, G., 2007. *Computational Science and Engineering*. Wellesley: Wellesley-Cambridge Press.

TRNSYS 16, 2007. *Mathematical Reference*. Madison, WI, USA: University of Wisconsin.

TRNSYS 17, 2009. *Mathematical reference*. Madison: Thermal Energy System Specialists, LLC.