Distributed Temperature Control via Geothermal Heat Pump Systems in Energy Efficient Buildings

TL;DR

This paper develops scalable, real-time distributed control algorithms for geothermal heat pump systems in multi-zone buildings, enhancing energy efficiency and user comfort without extensive measurements or predictions.

Contribution

It introduces a novel distributed control framework for multi-zone GHP systems that improves efficiency and scalability compared to traditional single-zone control schemes.

Findings

Control algorithms effectively balance comfort and energy use.

Distributed control does not require outdoor temperature predictions.

Simulation confirms improved efficiency and ease of implementation.

Abstract

Geothermal Heat Pump (GHP) systems are heating and cooling systems that use the ground as the temperature exchange medium. GHP systems are becoming more and more popular in recent years due to their high efficiency. Conventional control schemes of GHP systems are mainly designed for buildings with a single thermal zone. For large buildings with multiple thermal zones, those control schemes either lose efficiency or become costly to implement requiring a lot of real-time measurement, communication and computation. In this paper, we focus on developing energy efficient control schemes for GHP systems in buildings with multiple zones. We present a thermal dynamic model of a building equipped with a GHP system for floor heating/cooling and formulate the GHP system control problem as a resource allocation problem with the objective to maximize user comfort in different zones and to minimize the building energy consumption. We then propose real-time distributed algorithms to solve the control problem. Our distributed multi-zone control algorithms are scalable and do not need to measure or predict any exogenous disturbances such as the outdoor temperature and indoor heat gains. Thus, it is easy to implement them in practice. Simulation results demonstrate the effectiveness of the proposed control schemes.

Equations

$$\displaystyle(h(y))_{x}^{+}=\left\{\begin{array}[]{ll}h(y)&\textrm{if $x>0$}\\ \max(0,h(y))&\textrm{if $x=0$}\end{array}\right..$$

$$C_i{\dot{T}}_i=\frac{T^o-T_i}{R_i}+\sum _{j\in \mathcal{N}(i)}\frac{T_j-T_i}{R_{ij}}+\frac{T_{fi}-T_i}{R_{afi}}+Q_i$$

$$C_{fi}{\dot{T}}_{fi}=\frac{T_i-T_{fi}}{R_{afi}}+\frac{T_{wi}-T_{fi}}{R_{fwi}},i\in \mathcal{N}$$

$$C_{wi}{\dot{T}}_{wi}=\frac{T_{fi}-T_{wi}}{R_{fwi}}+c_w{q}_i(T_s-T_{fi}),i\in \mathcal{N}$$

$$\underset{Z_i,q_i,T_s,Z_{fi}}{\min }\sum _{i\in \mathcal{N}}\left[\frac{1}{2}{r}_i(Z_i-T_i^{set}{)}^2+s\frac{c_w{q}_i|T_s-Z_{fi}|}{-a{T}_s+b}\right]$$

$$\text{s. t. }\frac{T^o-Z_i}{R_i}+\sum _{j\in \mathcal{N}(i)}\frac{Z_j-Z_i}{R_{ij}}+\frac{Z_{fi}-Z_i}{R_{afi}}+Q_i=0$$

$$\frac{Z_i-Z_{fi}}{R_{afi}}+c_w{q}_i(T_s-Z_{fi})=0$$

$$0\le q_i\le q_i^{max}$$

$$T_s^{min}\le T_s\le T_s^{max}$$

$$C_i{\dot{T}}_i=\frac{T^o-T_i}{R_i}+\sum _{j\in \mathcal{N}(i)}\frac{T_{ij}-T_i}{R_{ij}}+\frac{T_{fi}-T_i}{R_{afi}}+Q_i$$

$$C_{ij}{\dot{T}}_{ij}=\frac{T_i-T_{ij}}{R_{ij}}+\frac{T_j-T_{ij}}{R_{ij}}$$

$$\frac{T^o-Z_i}{R_i}+\sum _{j\in \mathcal{N}(i)}\frac{Z_j-Z_i}{2R_{ij}}+\frac{Z_{fi}-Z_i}{R_{afi}}+Q_i=0$$

$$\underset{Z_i,u_i,Z_{fi}}{\min }\sum _{i\in \mathcal{N}}\left[\frac{1}{2}{r}_i(Z_i-T_i^{set}{)}^2+s\frac{|u_i|}{-a{T}_s+b}\right]$$

$$\text{s. t. }\frac{T^o-Z_i}{R_i}+\sum _{j\in \mathcal{N}(i)}\frac{Z_j-Z_i}{R_{ij}}+\frac{Z_{fi}-Z_i}{R_{afi}}+Q_i=0$$

$$\frac{Z_i-Z_{fi}}{R_{afi}}+u_i=0$$

$$0\le \frac{u_i}{c_w(T_s-Z_{fi})}\le q_i^{max}$$

$$L=$$

$$\displaystyle+\sum_{i\in\mathcal{N}}\zeta_{i}\Big{(}\frac{T^{o}-Z_{i}}{R_{i}}+\sum_{j\in\mathcal{N}(i)}\frac{Z_{j}-Z_{i}}{R_{ij}}+\frac{Z_{fi}-Z_{i}}{R_{afi}}+Q_{i}\Big{)}$$

$$\displaystyle+\sum_{i\in\mathcal{N}}\lambda_{i}\Big{(}\frac{Z_{i}-Z_{fi}}{R_{afi}}+u_{i}\Big{)}+\sum_{i\in\mathcal{N}}\mu_{i}^{-}(-u_{i})$$

$$+\sum _{i\in \mathcal{N}}{\mu }_i^{+}(u_i-q_i^{max}{c}_w(T_s-Z_{fi}))$$

$$\frac{\partial L}{\partial Z_i}=$$

$$+\sum _{j\in \mathcal{N}(i)}\frac{{\zeta }_j}{R_{ij}}+\frac{{\lambda }_i}{R_{afi}}=0,i\in \mathcal{N}$$

$$\frac{\partial L}{\partial u_i}=$$

$$\frac{\partial L}{\partial Z_{fi}}=$$

$$\frac{\partial L}{\partial {\zeta }_i}=$$

$$=$$

$$\frac{\partial L}{\partial {\lambda }_i}=$$

$${\mu }_i^{+}$$

$${\mu }_i^{-}$$

$${\dot{Z}}_i=$$

$$\displaystyle+\sum_{j\in\mathcal{N}(i)}\frac{1}{R_{ij}}\Big{)}-\sum_{j\in\mathcal{N}(i)}\frac{\zeta_{j}}{R_{ij}}-\frac{\lambda_{i}}{R_{afi}}\Big{)}$$

$${\dot{u}}_i=$$

$$\displaystyle-\mu_{i}^{+}+\mu_{i}^{-}+k_{eu_{i}}(\hat{u}_{i}-u_{i})\Big{)}$$

$${\dot{\hat{u}}}_i=$$

$${\dot{Z}}_{fi}=$$

$$\displaystyle-\mu_{i}^{+}q_{i}^{max}c_{w}+k_{eZ_{fi}}(\hat{Z}_{fi}-Z_{fi})\Big{)}$$

$${\dot{\hat{Z}}}_{fi}=$$

$${\dot{\zeta }}_i=$$

$$\displaystyle+\frac{Z_{fi}-Z_{i}}{R_{afi}}+Q_{i}\Big{)}$$

$${\dot{\lambda }}_i=$$

$${\dot{\mu }}_i^{+}=$$

$${\dot{\mu }}_i^{-}=$$

$$q_i=\frac{u_i}{c_w(T_s-Z_{fi})},i\in \mathcal{N}$$

$${\dot{\tilde{\zeta }}}_i=$$

$$+\frac{T_i-Z_i-T_{fi}+Z_{fi}}{R_{afi}},i\in \mathcal{N}.$$

$$\underset{Z_i,q_i,T_s,Z_{fi}}{\min }\sum _{i\in \mathcal{N}}\left[\frac{1}{2}{r}_i(Z_i-T_i^{set}{)}^2+s\frac{c_w^2{q}_i^2(T_s-Z_{fi}{)}^2}{-a{T}_s+b}\right]$$

Full text

I Introduction

According to an investigation by the United Nations, buildings are responsible for $40\%$ of energy consumption, $70\%$ of electricity consumption, and result in $30\%$ of greenhouse gas emission [1]. Roughly speaking, Heating Ventilation and Air Conditioning (HVAC) systems in buildings account for $40\%$ of the energy use [2]. It is therefore necessary to make them more energy efficient for environmental sustainability.

In recent years, Geothermal Heat Pump (GHP) systems are becoming popular among different HVAC systems, due to their highly efficient use of energy, i.e., they can usually deliver more than $3$kWh of heat with $1$kWh of electricity [3]. GHP systems are heating/cooling systems that use the ground as the temperature exchange medium. In winter, they transfer heat from the underground soil/water to buildings for heating, and vice versa in summer for cooling. Conventional control of the GHP system includes Proportional-Integral-Derivative (PID) control [4] and centralized Model Predictive Control (MPC) [5, 6]. These methods are practically efficient for cases with only a single thermal zone, however, they become either less efficient or costly (due to the centralized operation with heavy burdens of sensing, communication and computation of MPC) for cases with multiple thermal zones for large buildings. Since modern buildings are usually large, complex and are with multiple zones, scalable and easy-implementing control schemes are undoubtedly needed for them if equipped with GHP systems for heating/cooling.

This paper aims to develop real-time control schemes for GHP systems in typical multi-zone buildings. Specifically, we aim to design distributed algorithms to guide each controllable component to properly adapt their behavior such that system-wide objectives are achieved under given operating conditions. The emergence of distributed/decentralized control in network systems has been stimulated by smart sensing, communication, computing, and actuation technologies nowadays, e.g., in smart grids [7, 8], smart cities [9, 10], mobile robots [11], and intelligent transportation systems [12]. The advantages of distributed/decentralized control include: good scalability as the network grows; reduction of measurement, communication and computation compared with centralized control; privacy preserving. Thus, applying distributed/decentralized control to GHP system control and optimization is becoming an area of active research. Representative work includes, for example, distributed MPC [13, 3]. However, distributed MPC still requires a large amount of sensing, communication and computation. In most cases, it needs good prediction of future disturbances, i.e., outdoor temperature, sunlight, indoor occupancy, etc., which may be hard to obtain in reality. Different from the work using MPC, the controllers designed here are based on solving steady-state optimization problems via gradient algorithms: they (i) are dynamic feedback controllers that can be implemented without measuring or predicting disturbances, (ii) are scalable with respect to building structures, (iii) satisfy the system operating constraints, and (iv) ensure system efficiency, reliability and user comfort.

The structure of this paper is as follows. In Section II, we provide the detailed problem setup, including an introduction of the GHP system, a commonly used thermal dynamic model of the building network, and the optimization problem formulation. Since the original optimization problem is nonconvex, two different scenarios are considered in which the problem can be (approximated and) convexified: (i) the control inputs are only the water flow rates (Section III), and (ii) the control inputs are both water flow rates and the heat pump supply temperature (Section IV). For both scenarios, we use a modified primal-dual gradient method to design real-time distributed/decentralized control schemes. As a result, the thermal dynamics can be driven to equilibria which are the optimal solutions of those associated optimization problems. In Section V, two numerical examples are provided to illustrate the effectiveness of the designed control schemes, using a building with four adjacent zones. Finally, conclusions and future work are presented in Section VI.

Notation: The positive projection of a function $h(y)$ on a variable $x\in[0,+\infty)$, $(h(y))_{x}^{+}$ is:

[TABLE]

II Problem Setup

II-A GHP system in buildings

The schematic of a typical GHP system is illustrated in Figure 1, which consists of two hydronic and one refrigerant circuits, interconnected by two heat exchangers, i.e., an evaporator and a condenser [3, 6]. In the following we take the heating mode case as an example to explain the working process of the GHP system according to [14, 5], as the heating mode is more commonly used in practice.

The underground hydronic circuit contains a mixture of water and anti-freeze driven by a small circulating pump, and the temperature of its underground buried brine-filled side is relatively constant with a seasonal pattern, i.e., warmer in winter and cooler in summer than the outdoor temperature. The liquid refrigerant in the refrigerant circuit/the heat pump, first goes into the evaporator to absorb heat from the underground hydronic circuit and is converted to its gaseous state. Then this gaseous refrigerant passes through the compressor, and stops at the condenser in which it is converted to its liquid state to heat up the water in the distribution hydronic circuit. Finally this liquid refrigerant passes through the expansion valve, and stops at the evaporator for the next circulation. The distribution hydronic circuit is a grid of indoor under-surface pipes filled with water. Driven by another small circulating pump, this grid of pipes distributes heat to concrete floors (i.e., floor heating) or hydronic radiators (i.e., radiator heating) of a building for heating purpose. Here we only consider floor heating while radiator heating is similar and will be reported in a future paper.

In general, the heat pump consumes electrical power to transfer heat to the water in the distribution hydronic circuit. The amount of this heat depends on the flow rate of the water, the heat pump supply/forward temperature, and the return water temperature [15, 6]. Each zone in the building is equipped with a Thermal Wax Actuator (TWA) that adjusts the valve opening of the pipes for regulating the flow rate. The supply temperature is adjusted by regulating the compressor of the heat pump. The return water temperature can be approximated by the floor temperature which is accurate enough for control design [14]. These facts will be used later in system modeling as well as control design.

II-B Thermal dynamic model with a GHP system

According to the above configuration, we model a given building as an undirected connected graph $(\mathcal{N},\mathcal{E})$. Here $\mathcal{N}$ is the set of nodes representing zones/rooms, and $\mathcal{E}\subseteq\mathcal{N}\times\mathcal{N}$ is the set of edges. An edge $(i,j)\in\mathcal{E}$ means that zones $i$ and $j$ are neighbors. Let $\mathcal{N}(i)$ denote the set of neighboring zones of zone $i$. The thermal dynamics for each zone is described by a reduced Resistance-Capacitance (RC) model [16] (more discussion on this model is available in Remark 1):

[TABLE]

where $i\in\mathcal{N}$, $C_{i}$ is the thermal capacitance, $T_{i}$ is the indoor temperature, $T^{o}$ is the outdoor temperature, $R_{i}$ is the thermal resistance of the wall and window separating zone $i$ and outside, $R_{ij}$ is the thermal resistance of the wall separating zones $i$ and $j$, $T_{fi}$ is the floor temperature, $R_{afi}$ is the thermal resistance between the indoor air and the floor, and $Q_{i}\geq 0$ is the heat gain/disturbance from exogenous sources (e.g., user activity, solar radiation and device operation).

The thermal dynamics of floors equipped with a GHP system is described by a simplified lumped element model [6]:

[TABLE]

where $C_{fi}$ is the thermal capacitance of the floor, $T_{wi}$ is the temperature of the water in pipes, $R_{fwi}$ is the thermal resistance between the floor and the water, $C_{wi}$ is the thermal capacitance of the water, $c_{w}$ is the specific heat of the water, $q_{i}$ is the flow rate of the water, and $T_{s}$ is the supply temperature of the heat pump. Note that (i) $T_{s}$ is a common variable of the whole building [6], and (ii) the term $c_{w}q_{i}(T_{s}-T_{fi})$ stands for the heat transfer from the heat pump to the water in undersurface pipes [15].

Proposition 1**.**

When the GHP combined with floor heating/cooling system is off ($q_{i}=0$), (4)-(5) asymptotically converges to an equilibrium point which is uniquely determined by disturbances $T^{o},Q_{i}$. When the GHP combined with floor heating/cooling system is on, the asymptotic convergence property of (4)-(5) remains and the steady state is uniquely determined by disturbances $T^{o},Q_{i}$ and control inputs $q_{i},T_{s}$.

The above proposition can be directly derived by rearranging (4)-(5) in state-space representation, and showing that the system matrix is Hurwitz (an alternative way is to construct a quadratic Lyapunov function). So the desiderata is to design $q_{i},T_{s}$ only for periods when the GHP combined with floor heating/cooling system is on, more specifically, is to design the dynamics of $q_{i},T_{s}$ to drive (4)-(5) to some desired state.

II-C The optimization problem

In reality, each zone has a desired temperature which is the set point determined by users. The control objective considered in this paper is to regulate the temperature to be close to the set point in each zone, and to minimize the total energy consumption of the GHP system. More specifically, we consider the following steady-state optimization problem:

[TABLE]

where $i\in\mathcal{N}$ in (6b)-(6d), $r_{i}$ and $s$ are positive weight coefficients, $T_{i}^{set}$ is the temperature set point, $a,b$ are positive coefficients in the Coefficient of Performance (COP, i.e., an indicator of the relationship between the produced heat and consumed electricity by the heat pump [15]), $[0,q_{i}^{max}]$ is the range of $q_{i}$, and $[T_{s}^{min},T_{s}^{max}]$ is the range of $T_{s}$. Note that (i) to avoid confusion between steady-state values and temperature dynamics, we use $Z_{i},Z_{fi}$ to denote steady-state temperature values whereas $T_{i},T_{fi}$ are temperatures in the dynamic model (4)-(5), (ii) $T^{o}$ and $Q_{i}$ are exogenous disturbances, and (iii) we have merged the steady-state equations from (5) to obtain (6c). We assume that problem (6) is feasible and satisfies Slater’s condition [17]. Moreover, we have four important remarks.

$\bullet$ In the objective function (6a), the term relating to the total energy/electricity consumption (weighted by $s$) is given by $\sum_{i\in\mathcal{N}}\frac{c_{w}q_{i}|T_{s}-Z_{fi}|}{-aT_{s}+b}$: the term $\sum_{i\in\mathcal{N}}c_{w}q_{i}|T_{s}-Z_{fi}|$ stands for the total heat exchange between the heat pump and the water in pipes; $-aT_{s}+b>0$ is the COP which has been approximated as a linear function of $T_{s}$ as in [15]; $a,b$ can usually be obtained from the heat pump data sheet [15, 6].

$\bullet$ Parameters $r_{i},s$ are determined by users. If users prefer more comfort, they can increase $r_{i}$ and decrease $s$, and vice versa. Because of this flexibility, we do not impose constraints on the temperature comfortable range.

$\bullet$ In the heating mode, $T_{s}>Z_{fi}(\text{or }T_{fi}),\forall i$ hold; in the cooling mode, $T_{s}<Z_{fi}(\text{or }T_{fi}),\forall i$ hold. This is usually true in practice. For example [15, 6], in the heating mode, $T_{s}>27^{\circ}C$ while $T_{fi}<26^{\circ}C$. Once the mode is determined, the sign of $T_{s}-Z_{fi}(\text{or }T_{s}-T_{fi})$ is determined.

$\bullet$ The inequality constraints (6d)-(6e) are in accord with those in the optimization problem (7) in [6].

To conclude, the goal is to design the regulating rule for $q_{i},T_{s}$ so that system (4)-(5) can be driven to an equilibrium point which is the optimal solution to problem (6).

Remark 1**.**

[18]** Though system (4) is a 1st-order RC model, using higher order RC models does not affect the formulation of (6) since it is a steady-state optimization problem. For example, for the 2nd-order model in [19, 16]

[TABLE]

where $T_{ij}$ is the temperature of the wall separating zones $i$ and $j$, and $C_{ij}$ is the thermal capacitance of the wall, the corresponding steady-state Equation (6b) is given by

[TABLE]

which results from the steady-state equation $Z_{ij}=\frac{Z_{i}+Z_{j}}{2}$ ($Z_{ij}$ is the steady-state temperature of the wall), i.e., the steady-state equation of the higher order model can be reduced to (6b) by eliminating states of solids in the building envelope. Since our control design procedures proposed later are based on solving the steady-state optimization problem (6), using higher order models will not affect them.

III Scenario I: Flow Rate Control Only

III-A Problem reformulation

In this section, we consider the water flow rate $q_{i}$ as the only control input to each zone, and regard $T_{s}$ as a known exogenous signal which satisfies constraint (6e). Such scenario could happen in practice, for instance, $T_{s}$ is required to track some prescribed curve [3]. In this case, problem (6) can be simplified into

[TABLE]

where $i\in\mathcal{N}$ in (7b)-(7d), $u_{i}=c_{w}q_{i}(T_{s}-Z_{fi})$ is introduced to replace terms on $q_{i}$, and constraint (6e) is dropped (under these actions, problems (6) and (7) are still equivalent). Note that the GHP system is in either the heating mode or the cooling mode. Once the mode is determined, the signs of $u_{i}$ and $T_{s}-Z_{fi}$ are determined so that (i) $|u_{i}|$ equals either $u_{i}$ or $-u_{i}$, and (ii) the inequality constraint (7d) can become linear by multiplying $c_{w}(T_{s}-Z_{fi})$ on both sides. Thus, problem (7) naturally becomes convex.

III-B A distributed algorithm

Once the mode is determined, problem (7) can be solved in either a centralized or distributed/decentralized way. Any centralized algorithm requires to measure the outdoor temperature $T^{o}$ and the indoor heat gain $Q_{i}$ in every zone ($T_{s}$ is given). Because these exogenous disturbances can fluctuate frequently and are not easy to obtain, the cost of centralized algorithms would be expensive. Next we develop a real-time distributed algorithm that does not need measurement of these exogenous disturbances.

Consider the heating mode case (the cooling mode case is similar). The Lagrangian function of (7) is given by

[TABLE]

where $\zeta_{i},\lambda_{i},\mu_{i}^{+},\mu_{i}^{-}$ are the Lagrange multipliers/dual variables for constraints (7b)-(7d). Since problem (7) is convex, feasible and satisfies Slater’s condition, the Karush-Kuhn-Tucker (KKT) conditions are necessary and sufficient conditions for optimality [17], given by

[TABLE]

Motivated by a modified primal-dual gradient method [20, 21], we design the following algorithm to solve (7):

[TABLE]

where $i\in\mathcal{N}$, $k_{Z_{i}},k_{u_{i}},k_{eu_{i}},\hat{k}_{eu_{i}},k_{Z_{fi}},k_{eZ_{fi}},\hat{k}_{eZ_{fi}},k_{\zeta_{i}},k_{\lambda_{i}},$ $k_{\mu_{i}^{+}},k_{\mu_{i}^{-}}$ are positive scalars representing the controller gains, and we have introduced the auxiliary states $\hat{u}_{i},\hat{Z}_{fi}$: since the objective function (7a) is not strictly convex in $u_{i},Z_{fi}$, a standard primal-dual gradient method [22] could yield large oscillations; after introducing the extra dynamics, the transient behavior of the overall system can be improved (demonstrated in Section V). Note that $T_{i},T_{fi}$ have their own dynamics given by (4)-(5) and thus can not be designed, which is why we replace $T_{i},T_{fi}$ with $Z_{i},Z_{fi}$ initially, i.e., $Z_{i},Z_{fi},i\in\mathcal{N}$ are ancillary state variables. According to [21, 23], it is true that (9) asymptotically converges to an equilibrium point which is the optimal solution of (7), since the optimization problem is convex and extra dynamics have been included in (9). Now using

[TABLE]

as the control input to system (4)-(5), we can naturally obtain a real-time distributed controller to regulate (4)-(5) to a steady state which is the optimal solution to problem (7) ((7) is equivalent to (6) under known $T_{s}$).

Theorem 1**.**

Given constant/step change/slow-varying $T^{o},$ $Q_{i},T_{s}$, the trajectory of system (4)-(5) and (9)-(10) asymptotically converges to an equilibrium point at which $T_{i},q_{i},T_{fi}$ of the equilibrium point is the optimal solution of (6).

Proof.

According to [21, 23], each trajectory of system (9) asymptotically converges to an equilibrium point which is the optimal solution of (7). Under (10), the resulting equilibrium point after state transformation from $u_{i}$ to $q_{i}$, is the optimal solution of (6). On the other hand, the trajectory of the overall system (4)-(5) and (9)-(10) also asymptotically converges to an equilibrium point, due to the cascade nature, i.e., (9)-(10)$\rightarrow$(4)-(5). By Proposition 1, since the equilibrium point of (4)-(5) is uniquely determined by the inputs $T^{o},Q_{i},q_{i},T_{s}$ in which $q_{i}$ is given by (10), we have $T_{i}=Z_{i},T_{fi}=Z_{fi}$ (here $T_{i},T_{fi}$ are states given by (4)-(5)) when the overall system reaches steady state, which completes the proof. ∎

This theorem requires $T^{o},Q_{i}$ to be either constant, step change, or slow-varying, which holds in practice as they vary at a time-scale of minutes. Remark that our controller operates in real-time, i.e., at a time-scale of seconds.

In Equation (9f), the disturbances $T^{o},Q_{i}$ appear. Motivated by [21], to make the algorithm implementable without measuring these terms, we introduce $\tilde{\zeta}_{i}=\frac{\zeta_{i}}{k_{\zeta_{i}}}-C_{i}T_{i}$ as

[TABLE]

Moreover, we substitute $\zeta_{i}=k_{\zeta_{i}}(\tilde{\zeta}_{i}+C_{i}T_{i})$ into (9a), (9d) to eliminate $\zeta_{i}$. Now the proposed control scheme (9a)-(9e), (9g)-(9i) and (10)-(III-B) is completely distributed as shown in Figure 2 and can be implemented as follow. Given $T_{s},C_{i},R_{i},R_{ij},R_{afi},r_{i},s,a,b,q_{i}^{max}$, each zone in the building collects $T_{i}^{set}$ from users, locally measures its indoor temperature $T_{i}$ and floor temperature $T_{fi}$, receives the feedback signals $\zeta_{j}=k_{\zeta_{j}}(\tilde{\zeta}_{j}+C_{j}T_{j})$ and $T_{j}-Z_{j}$ from its neighboring zones, and then uses the information to update $Z_{i},u_{i},\hat{u}_{i},Z_{fi},\hat{Z}_{fi},\tilde{\zeta}_{i},\lambda_{i},\mu_{i}^{+},\mu_{i}^{-},q_{i}$. Here $C_{i},R_{i},R_{ij},R_{afi},a,b,q_{i}^{max}$ are building parameters, $r_{i},s$ are parameters specified by users, and $T_{s}$ is a known signal.

Remark 2**.**

In reality, due to that $Z_{i},Z_{j}$ (or $T_{i},T_{j}$) of neighboring zones are often very close to each other and $R_{ij}$ is not small, the total heat gain/loss from neighboring zones is (sometimes much) less dominant compared with the heat gain/loss from the outside plus the indoor heat gain in every zone. Thus, the term $\sum_{j\in\mathcal{N}(i)}\frac{T_{j}-T_{i}}{R_{ij}}$ could be ignored in (4) as well as the term $\sum_{j\in\mathcal{N}(i)}\frac{Z_{j}-Z_{i}}{R_{ij}}$ in (6)/(7). In this case, following the same design procedure, we end up with a completely decentralized control scheme (given by setting $R_{ij}=\infty$ in (9a) and (III-B)) that only needs local measurement.

IV Scenario II: The General Case

IV-A Problem reformulation

In this section, both $q_{i}$ and $T_{s}$ are considered as control inputs to the system. Instead of handling the nonconvex problem (6) directly, we focus on an approximate version of this problem, i.e., rather than minimizing the exact energy consumption in the objective function, we minimize one of its upper bound. This approximation can help convexify (6), which will become clear later. The approximate problem is

[TABLE]

where we have modified the term $\sum_{i\in\mathcal{N}}\frac{c_{w}q_{i}|T_{s}-Z_{fi}|}{-aT_{s}+b}$ in (6a) to be $\sum_{i\in\mathcal{N}}\frac{c_{w}^{2}q_{i}^{2}(T_{s}-Z_{fi})^{2}}{-aT_{s}+b}$ in (12a). The latter term is actually an upper bound of the square of the exact energy consumption: using a Cauchy-Schwarz inequality, we have

[TABLE]

where $|\mathcal{N}|$ is the number of zones in the building. Therefore, rather than minimizing the total consumption directly, the objective here aims to minimize its upper bound, which is sufficient for energy saving purpose.

Now we show how (12) can be turned convex. Again, we introduce variables $u_{i}=c_{w}q_{i}(T_{s}-Z_{fi}),i\in\mathcal{N}$ to get

[TABLE]

Note that (i) the function $E=\frac{\sum_{i\in\mathcal{N}}u_{i}^{2}}{-aT_{s}+b}$ is convex in $u_{i},T_{s}$ as its Hessian matrix equals

[TABLE]

which is positive semi-definite and (ii) once the mode is determined, the signs of $u_{i}$ and $T_{s}-Z_{fi}$ are determined so that constraint (7d) can become linear by multiplying $c_{w}(T_{s}-Z_{fi})$ on both sides. Next we design a distributed algorithm to solve the convex optimization problem (13).

IV-B A distributed algorithm

Similar to Section III-B, we consider the heating mode case. The design methodology is motivated by a modified primal-dual gradient method [20, 21]. For simplicity, we directly write down the resulting distributed algorithm:

[TABLE]

where $i\in\mathcal{N}$, $k_{Z_{i}},k_{u_{i}},k_{T_{s}},k_{Z_{fi}},k_{eZ_{fi}},\hat{k}_{eZ_{fi}},k_{\zeta_{i}},k_{\lambda_{i}},k_{\mu_{i}^{+}},$ $k_{\mu_{i}^{-}},k_{\nu^{+}},k_{\nu^{-}}$ are positive scalars representing the controller gains, and we have introduced the auxiliary state $\hat{Z}_{fi}$ to improve the performance of the algorithm (since $u_{i}$ and $T_{s}$ are coupled in $E$ whose Hessian matrix is not always zero, only adding $\hat{Z}_{fi}$ is enough for behavior enhancement). Now using (10) and (19c) as the control input to (4)-(5) where $u_{i},T_{s},Z_{fi}$ are determined by (19), we can naturally obtain a real-time distributed controller to regulate (4)-(5) to a steady state which is the optimal solution to (12).

Theorem 2**.**

Given constant/step change/slow-varying $T^{o},$ $Q_{i}$ (remark that they vary at a time-scale of minutes), each trajectory of the overall system (4)-(5), (10) and (19) asymptotically converges to an equilibrium point at which $T_{i},q_{i},T_{s},T_{fi}$ of this point is the optimal solution of (12).

Proof.

The proof is similar to that of Theorem 1, and thus, is omitted for brevity. ∎

In Equation (19f), the disturbances $T^{o},Q_{i}$ appear. Similar to Section III-B, to make the algorithm implementable without measuring these terms, we introduce $\tilde{\zeta}_{i}=\frac{\zeta_{i}}{k_{\zeta_{i}}}-C_{i}T_{i}$ whose dynamics is given by (III-B). Moreover, we substitute $\zeta_{i}=k_{\zeta_{i}}(\tilde{\zeta}_{i}+C_{i}T_{i})$ into (19a), (19d) to eliminate $\zeta_{i}$. Now the proposed control algorithm (19a)-(19e), (19g)-(19k) and (10)-(III-B) is completely distributed as shown in Figure 3 and can be implemented as follow. Given $C_{i},R_{i},R_{ij},R_{afi},r_{i},s,a,b,q_{i}^{max}$, each zone in the building collects $T_{i}^{set}$ from users, locally measures its indoor temperature $T_{i}$ and floor temperature $T_{fi}$, receives the feedback signals $\zeta_{j}=k_{\zeta_{j}}(\tilde{\zeta}_{j}+C_{j}T_{j}),T_{j}-Z_{j}$ from its neighboring zones and $T_{s}$ from the compressor, and then uses the information to update $Z_{i},u_{i},Z_{fi},\hat{Z}_{fi},\tilde{\zeta}_{i},\lambda_{i},\mu_{i}^{+},\mu_{i}^{-},q_{i}$. On the other hand, given $T_{s}^{min},T_{s}^{max}$, the compressor receives the feedback signals $u_{i},\mu_{i}^{+}q_{i}^{max}c_{w}$ from each zone, updates $T_{s},\nu^{+},\nu^{-}$, and then broadcasts $T_{s}$. Here $C_{i},R_{i},R_{ij},R_{afi},a,b,q_{i}^{max},T_{s}^{min},T_{s}^{max}$ are building parameters, and $r_{i},s$ are specified by users.

Remark 3**.**

Similar to Remark 2, the term $\sum_{j\in\mathcal{N}(i)}\frac{T_{j}-T_{i}}{R_{ij}}$ could be ignored in (4) as well as the term $\sum_{j\in\mathcal{N}(i)}\frac{Z_{j}-Z_{i}}{R_{ij}}$ in (12). Then following the same design procedure, we obtain another distributed control scheme (given by setting $R_{ij}=\infty$ in (19a) and (III-B)) that does not need communication between neighboring zones, i.e., it requires less communication.

V Numerical Investigations

In this section, we present two numerical examples for scenarios described in Sections III and IV respectively, using a house with four adjacent zones as illustrated in Figure 1. Only the heating case is presented in the following, while the cooling case is similar under the proposed control schemes.

The parameters of the simulations are obtained from [15, 24, 25]: all $C_{i}=20$kJ$\text{/}^{\circ}$C, all $C_{fi}=35$kJ$\text{/}^{\circ}$C, all $C_{wi}=25$kJ$\text{/}^{\circ}$C, all $R_{i}=15^{\circ}$C/kW, all $R_{ij}=23^{\circ}$C/kW, all $R_{afi}=3^{\circ}$C/kW, all $R_{fwi}=5^{\circ}$C/kW, $c_{w}=4.186$kJ/kg$\text{/}^{\circ}$C, $[q_{i}^{max}]=[0.03,0.04,0.045,0.035]$kg/s, $[T_{s}^{min},T_{s}^{max}]=[38,42]^{\circ}$C, $a=0.11\text{/}^{\circ}$C, $b=8.4$, all $r_{i}=0.5$p.u., all $k_{Z_{i}}=0.025$p.u., all $k_{Z_{fi}}=0.033$p.u., $k_{T_{s}}=0.05$p.u., all $k_{u_{i}}=\hat{k}_{eZ_{fi}}=k_{\zeta_{i}}=k_{\lambda_{i}}=k_{\mu_{i}^{+}}=k_{\mu_{i}^{-}}=k_{\nu^{+}}=k_{\nu^{-}}=1$p.u., all $k_{eu_{i}}=10$p.u., all $\hat{k}_{eu_{i}}=0.1$p.u., and all $k_{eZ_{fi}}=2$p.u. (p.u. means per unit). The outdoor temperature, indoor heat gains, and supply temperature in Scenario I are shown in Figure 4.

The simulation result of the first scenario is illustrated in Figures 5-6, in which we set $s=0$ before $15$h and $s=10$p.u. thereafter. The curves labelled with “app” indicate the case of using the decentralized controller given in Remark 2, i.e., a communication free scheme. It can be seen that the difference between these two cases is not large, i.e., less than $1.3^{\circ}$C in temperatures, indicating that the performance of the decentralized controller could be acceptable in practice. Before $15$h, since there is no consumption reduction purpose, i.e., $s=0$, the temperature trajectories under (9a)-(9e), (9g)-(9i) and (10)-(III-B) track their set points unless the corresponding water flow rate saturates (although not shown here, using a PID controller will result in the same temperature trajectories during this period). After $15$h, deviations from temperature set points appear due to the consideration of energy saving (while only using a PID controller can not reduce energy consumption unless forcing users to change their set points). On the other hand, we compare the flow rate response under the distributed controller, with and without extra dynamics. The result shown in Figure 7 demonstrates that after introducing those extra dynamics, the system performance has been improved that the oscillations are largely attenuated.

In the second scenario, $s=1$p.u. holds before $13$h and $s=5$p.u. thereafter. We can see that the difference between using the distributed controller (19a)-(19e), (19g)-(19k) and (10)-(III-B) and using its simplified version given in Remark 3 is also not large, as shown in Figures 8-10. The temperature deviations with respect to their set points before increasing the weight coefficient $s$ are smaller than those thereafter since starting from $13$h, energy saving becomes more important while user comfort becomes less. All these two scenarios inspire us that tuning the weight coefficient $s$ (or equivalently $r_{i}$, as the optimal solution of (6)/(12) depends on the ratio $r_{i}/s,i\in\mathcal{N}$) can balance user comfort and GHP system energy consumption – there always exists a tradeoff between user comfort and energy saving.

VI Conclusion and Future Work

This paper presents distributed control frameworks on real-time temperature regulation via GHP combined with floor heating/cooling systems in energy efficient buildings. The controllers regulate water flow rates as well as the heat pump supply temperature, which balance user comfort and energy saving. Moreover, they can automatically adapt to changes of disturbances such as the outdoor temperature and indoor heat gains, without measuring or predicting those values. Also, the implementation of the controllers are simple.

Future work includes: developing distributed/decentralized control schemes to drive system (4)-(5) to the exact optimal solution of problem (6) (this may require an exact convex relaxation approach to solve (6)); considering both floor and radiator heating/cooling (this may require model reduction for the dynamics of hydronic radiators); and investigating, e.g., the $H_{2}$ and $H_{\infty}$ performances of the controlled systems.