An integrated model predictive control approach for optimal HVAC and energy storage operation in large-scale buildings

TL;DR

This paper presents a specialized Model Predictive Control strategy for cost-effective operation of large-scale smart buildings, integrating HVAC, energy storage, photovoltaic generation, and demand-response participation.

Contribution

It introduces a computationally efficient MPC approach for managing HVAC and storage in large-scale buildings considering multiple energy sources and demand-response constraints.

Findings

The proposed MPC achieves cost savings in building operation.

The approach remains computationally feasible for large-scale systems.

Performance is validated under uncertainties and disturbances.

Abstract

This paper deals with the problem of cost-optimal operation of smart buildings that integrate a centralized HVAC system, photovoltaic generation and both thermal and electrical storage devices. Building participation in a Demand-Response program is also considered. The proposed solution is based on a specialized Model Predictive Control strategy to optimally manage the HVAC system and the storage devices under thermal comfort and technological constraints. The related optimization problems turn out to be computationally appealing, even for large-scale problem instances. Performance evaluation, also in the presence of uncertainties and disturbances, is carried out using a realistic simulation framework.

Tables

Table 1. Table 1: Test building characteristics

Building Component		Value
Weather and Location		Turin (Italy)
Floor Area [$m^2$]		375
Floor [#]		8
Zone [#]		126
Window to Wall Ratio [%]		20
Solar Transmittance [%]		30
Internal Loads	Occupants [#]	120
	Lighting [$W/m^2$]	4.00
	Equipment [$W/m^2$]	3.25
Heating: Heat Pump	Heating Capacity [$kW$]	389
	Heating Power Consumption [$kW$]	82
	TES Volume [$m^3$]	8
Cooling: Heat Pump	Cooling Capacity [$kW$]	389
	Cooling Power Consumption [$kW$]	82
Electric Energy Storage	Capacity [$kWh$]	28
	Max Charging Rate [$kW$]	24
	Max Discharging Rate [$kW$]	24
PV Plant	Panels [#]	5
	Plant Surface Area [$m^2$]	100
	Plant Peak Power [$kW$]	10

Table 2. Table 2: Test building construction materials (name/thickness [mm])

	External Walls	Internal Walls
Outside Layer	Brickwork/100	Gypsum plaster/13
Layer 2	Extruded polystyrene/80	Brickwork/10
Layer 3	Concrete block/100	Gypsum plaster/13
Layer 4	Gypsum plaster/15

Table 3. Table 3: Winter season - DR program

	$h_j$	${\mu }_j$	$S_j$ [$kWh$]	$R_j$ [€]
${ℛ}_1$	107	5	37	3.6
${ℛ}_2$	178	6	40	4.4
$ℛ3$	325	3	6	2.8

Table 4. Table 4: Heating operation mode - Simulation results over 3 days

Cost without DR [€]	136.14
No. of fulfilled DR requests	2
DR reward [€]	8.00
Overall cost with DR [€]	128.14
Worst zone average bound violation [°C]	0.139

Table 5. Table 5: Heating operation mode simulation - Comparison with thermostatic control (without DR)

	Proposed MPC	Thermostatic Control
Overall cost [€]	136.64	183.39
Worst zone average bound violation [°C]	0.139	0.142

Table 6. Table 6: Cooling operation mode - Simulation results over 3 days

Cost without DR [€]	79.68
No. of fulfilled DR requests	3
DR reward [€]	6.40
Overall cost with DR [€]	73.28
Worst zone average bound violation [°C]	0.145

Table 7. Table 7: Cooling operation mode simulation - Comparison with thermostatic control (without DR)

	Proposed MPC	Thermostatic Control
Overall cost [€]	76.74	119.29
Worst zone average bound violation [°C]	0.156	0.169

Table 8. Table 8: Simulation results for exact and inaccurate forecasts (worst-case over 20 realizations)

	Overall cost [€]	Worst zone average bound violation [°C]
Exact forecast	136.64	0.139
Uncertain forecast on external temperature	144.82	0.141
Uncertain forecast on solar irradiance	144.98	0.144
Uncertain forecast on zone internal gains	143.90	0.144
Uncertain forecast on all three inputs	147.28	0.146

Equations

$$\mathbf{T}(k+1)={\Theta }^{\mathbf{T}}{\mathbf{\Phi }}^{\mathbf{T}}(k),$$

$$\begin{array}[]{r}\mathbf{\Phi}^{\mathbf{T}}(k)=[\mathbf{T}^{\prime}(k)~{}\dots~{}\mathbf{T}^{\prime}(k-k_{\mathbf{T}})\\ \mathbf{h}^{\prime}(k)~{}\dots~{}\mathbf{h}^{\prime}(k-k_{\mathbf{h}})\\ \mathbf{e}^{\prime}(k)~{}\dots~{}\mathbf{e}^{\prime}(k-k_{\mathbf{e}})]^{\prime},\end{array}$$

$$\underline{\mathbf{T}}(k)\le \mathbf{T}(k)\le \overline{\mathbf{T}}(k).$$

$${\mathcal{C}}^T(k)=\{\text{mbox}\text{eqref}eq:Tdyn,\text{eqref}eq:Tdyn2\},\mathcal{C}(k)=\{\text{mbox}\text{eqref}eq:comfort\}.$$

$$0\le v_i(k)\le {\overline{v}}_i.$$

$$h_i(k)={\gamma }_i\left(T^{SND}(k)-T_i(k)\right)v_i(k),$$

$${\mathcal{C}}_H^h(k)=\{0\le h_i(k)\le {\gamma }_i\left(T^{SND}(k)-T_i(k)\right){\overline{v}}_i,i=1,\dots ,m\}$$

$${\mathcal{C}}_C^h(k)=\{0\ge h_i(k)\ge {\gamma }_i\left(T^{SND}(k)-T_i(k)\right){\overline{v}}_i,i=1,\dots ,m\}$$

$$T^{SND}(k)=T^{TES}(k).$$

$$T^{H{P}_H}(k+1)=T_0^{H{P}_H}(k),$$

$$T_{in}^{H{P}_H}(k)=T^{TES}(k).$$

$$T_0^{H{P}_H}(k)>T_{in}^{H{P}_H}(k)=T^{TES}(k),$$

$$W^{H{P}_H}(k)={\alpha }^{H{P}_H}\left(T_0^{H{P}_H}(k)-T^{TES}(k)\right),$$

$$T_0^{H{P}_H}(k)\ge T^{TES}(k)$$

$$T^{TES}(k+1)={\Theta }^{TES}{\mathbf{\Phi }}^{TES}(k),$$

$$\begin{array}[]{r}\mathbf{\Phi}^{TES}(k)=[T^{TES}(k)~{}\dots~{}T^{TES}(k-k_{TT})\\ T^{HP_{H}}(k)~{}\dots~{}T^{HP_{H}}(k-k_{TH})\\ \mathbf{h}^{\prime}(k)~{}\dots~{}\mathbf{h}^{\prime}(k-k_{T\mathbf{h}})]^{\prime},\end{array}$$

$$T^{H{P}_C}(k+1)=T_0^{H{P}_C}(k).$$

$$W^{H{P}_C}(k)={\alpha }^{H{P}_C}\left(T_{in}^{H{P}_C}(k)-T_0^{H{P}_C}(k)\right),$$

$$T_0^{H{P}_C}(k)\le T_{in}^{H{P}_C}(k).$$

$$T_{in}^{H{P}_C}(k+1)={\Theta }^C{\mathbf{\Phi }}^C(k),$$

$$\begin{array}[]{r}\mathbf{\Phi}^{C}(k)=[T^{HP_{C}}_{in}(k)~{}\dots~{}T^{HP_{C}}_{in}(k-k_{TC})\\ T^{HP_{C}}(k)~{}\dots~{}T^{HP_{C}}(k-k_{HC})\\ \mathbf{h}^{\prime}(k)~{}\dots~{}\mathbf{h}^{\prime}(k-k_{C\mathbf{h}})]^{\prime}\end{array}$$

$${\mathcal{C}}_H(k)=\{\text{mbox}\text{eqref}eq:HPH,\text{eqref}eq:HPpower,\text{eqref}eq:HPactivation,\text{eqref}eq:Hdynamics1,\text{eqref}eq:Hdynamics2\},{\mathcal{C}}_C(k)=\{\text{mbox}\text{eqref}eq:HPC,\text{eqref}eq:HPpowerC,\text{eqref}eq:HPactivationC,\text{eqref}eq:Cdynamics1,\text{eqref}eq:Cdynamics2\}.$$

$$E^{EES}(k+1)=E^{EES}(k)+\eta W_{+}^{EES}(k)-\frac{1}{\eta }{W}_{-}^{EES}(k),$$

$$W^{EES}(k)=W_{+}^{EES}(k)-W_{-}^{EES}(k).$$

$$0\le W_{+}^{EES}(k)\le {\overline{W}}_{+}^{EES},$$

$$0\le W_{-}^{EES}(k)\le {\overline{W}}_{-}^{EES},$$

$$0\le E^{EES}(k)\le {\overline{E}}^{EES},$$

$${\mathcal{C}}^{EES}(k)=\{\text{mbox}\text{eqref}eq:EESdyn,\text{eqref}eq:EESdyn2,\text{eqref}eq:EESconstr1,\text{eqref}eq:EESconstr11,\text{eqref}eq:EESconstr2\}.$$

$${\overline{W}}^{PV}(k)={\theta }_1^{PV}I(k)+{\theta }_2^{PV}{I}^2(k)+{\theta }_3^{PV}I(k)T^A(k),$$

$${\overline{W}}^{PV}(k)={\Theta }^{PV}\mathbf{e}(k)$$

Full text

††journal: Applied Energy

1 Introduction

It is well-known that buildings represent one of the major electricity consumers worldwide. About a half of the total energy consumption is devoted to heating, ventilation and air conditioning (HVAC). While HVAC devices in older buildings are operated through simple rules, newer buildings exploit the opportunities offered by information and communication technologies to efficiently operate such appliances [1]. An efficient control of most of the building devices yields to a global reduction of energy consumption, with a consequent reduction of emissions and energy bills [2].

Although a large part of the currently deployed building management systems are rule-based, Model Predictive Control (MPC) is gaining a lot of importance, owing to its flexibility and its ability to take a number of different requirements and constraints into account [3]. Indeed, to optimize the building operation cost, several applications of MPC can be found in the literature, where both linear [4] and non-linear [5] dynamics are considered. In [6], an MPC-based enthalpy control algorithm has been derived. Other works deal with uncertainty in models and/or exogenous variables by using stochastic MPC approaches. Specifically, in [7, 8] a tractable reformulation using chance constraints has been devised, while a scenario-based approach has been exploited in [9].

In order to reduce the computational complexity which often affects MPC implementations, in [10] a machine learning approach based on multivariate regression and dimensionality reduction algorithms has been proposed, and a case study involving a 6-zone building has been carried out. In [11], a cooperative Fuzzy MPC has been implemented on a real building composed of five zones.

Due to the relevance of buildings as players in smart grids, an important aspect recently considered in the building management literature is participation in Demand-Response (DR) programs, which enables end-users to become active players in the electricity system [12, 13, 14]. DR is usually managed by an intermediary player (known as the aggregator [15]) whose role is to gather flexibility from its affiliated consumers, in order to build services to be sold either to the grid operator or on the wholesale market. In response to the aggregator’s requests, which become known in advance, consumers may choose to adjust their consumption patterns in order to comply and receive a monetary reward [16, 17]. Reduction of energy cost in presence of DR has been studied in [18] via an agent-based modeling approach.

Considering the growing importance of distributed generation and the relevant share of energy consumption by buildings, energy storage systems are expected to become increasingly used. Such systems are employed to achieve high levels of self-sufficiency, as well as system resiliency [19, 20]. Furthermore, optimal sizing and efficient management of electrical storage in smart buildings, mainly in the presence of photovoltaic (PV) generation [21] and plug-in electric vehicles [22], is a widely investigated topic. In view of the participation of a building in DR programs, the presence of energy storage facilities contributes to achieve a high level of consumption flexibility.

1.1 Paper contribution

The aim of this work is to propose a novel MPC-based control strategy aimed at the minimization of the electricity bill in a large-scale building. The building is assumed to be equipped with a central HVAC system powered by heat pumps (HPs), a PV plant, thermal (TES) and electrical (EES) storage devices, and to participate in a DR program. In the proposed approach, heat pumps, heater actuators, storage devices, and PV generation facilities are managed in an optimal fashion by a centralized control unit. The employed DR paradigm consists of price-volume signals sent by an aggregator, which specify a monetary reward granted to the building operator upon the fulfillment of energy consumption constraints within a given time period [23, 24].

The control law for the different actuating devices of the overall HVAC system is computed through a receding horizon algorithm involving a two-step optimization strategy. Several constraints are taken into account, namely thermal comfort, DR programs, operating limits of heat exchangers, PV power availability, storage capacity and maximum charging/discharging rates. Meteorological variables, building occupancy, light and internal appliance loads are considered as exogenous disturbances, of which the availability of measured/predicted values is assumed.

The main contributions of the paper in the above context are as follows.

• •

Since the computational complexity of the control algorithm is a crucial aspect when dealing with large-scale buildings, a framework is adopted in which heat flows in the different building zones are regulated by controlling the air mass flow through each fan coil. This choice, along with approximate linear modeling of the building zones, allows for formulating the optimal control problem at each step as a Linear Program (LP). Hence, the involved computational burden scales nicely with the problem dimension and the solver provides a fast solution even for a large number of optimization variables and constraints [25]. Actually, the presence of DR programs calls for the introduction of binary decision variables, changing the nature of the problem from LP to Mixed Integer Linear Program (MILP). Nevertheless, the number of such variables amounts at most to few units, i.e., the number of time intervals that contain a DR request falling inside the optimization horizon and, more importantly, it is independent of the building size (i.e., the number of zones). As a result, the computational burden of the overall optimization algorithm is affordable even for large-scale instances, making the proposed approach suitable for real-world applications with hundreds of zones.

• •

As it is well-known in the literature, multiple sources of uncertainty come into the picture when dealing with building heating/cooling control, such as modeling inaccuracies, disturbances, weather forecasting errors and possibly undersized HVAC system. For this reason, when implementing a constrained MPC strategy, the presence of such uncertainties often causes infeasibility of the optimization problem. A standard way to deal with this issue is to replace hard bounds with soft bounds on the constraints and introduce penalties in the cost function accounting for constraint violations. Unfortunately, the choice of appropriate weights for the different penalty terms becomes a formidable task, especially in view of the high number of variables/constraints involved and the hard-to-predict “size” of constraint violations. To overcome these issues, this paper proposes a different approach based on a two-step optimization strategy. At each time step, first an ancillary LP is solved in order to reset the comfort constraints, by minimizing the $l_{1}$ norm of the bound violations. Then, the feasible overall optimization problem is solved using the modified constraints. A further advantage of this approach is that it allows for quick recovery from realistic situations in which the original comfort constraints cannot be satisfied.

To validate the proposed control algorithm, numerical simulations have been performed on EnergyPlus [26], an industry-standard realistic building model simulator. In particular, performance analysis in the presence of uncertain forecasts of exogenous variables has been carried out.

The paper is organized as follows: in Section 2 an overview of the proposed control architecture is presented, in Section 3 the relevant models are introduced. In Section 4 the control problem is formulated, its solution derived and implementation issues discussed. The test cases are presented in Section 5 and discussed in Section 6. Finally, conclusions are drawn in Section 7.

[TABLE]

[TABLE]

[TABLE]

[TABLE]

2 System architecture overview

This paper focuses on a building composed of several zones and equipped with a central HVAC system. Heating and cooling power is provided by electrical heat pumps. The heating pump HPH is assumed to be connected to a thermal storage device (TES), while the cooling pump HPC is not. Although thermal storage is quite common for both heating and cooling use, it is here employed in the first case only, with the aim of illustrating the performance of the proposed control strategy both in the presence and in the absence of a storage device. Heat exchange at zone level is provided by fan coil units, one for each zone, exchanging heat with the TES at constant water throughput. Heat flow regulation is achieved by controlling the air mass flow through each fan coil. It is also assumed that each zone is equipped with a temperature sensor. The building also features short-term electrical energy storage (EES) facilities (e.g., via a rechargeable battery system) and PV generation. Participation of the building in a Demand-Response (DR) program is also assumed. Regulation of the HPs, TES, fan coil air flow, EES and PV is implemented via a digital centralized control unit that operates in discrete time with sampling period $\tau_{s}$.

3 System modeling

A block diagram of the overall control system is depicted in Figure 1 along with the relevant signals. The modeling of each component is detailed in the following subsections.

3.1 Zone thermal model

The $m$ building zones are denoted by ${\rm Z}_{1},\dots,{\rm Z}_{m}$, and the respective heat exchanger devices by ${\rm H}_{1},\dots,{\rm H}_{m}$. All zones are equipped with temperature sensors, and the air flow of the heat exchangers can be independently regulated by the control unit. The control action must ensure that each room temperature satisfies a time-varying comfort constraint.

Let us define the following variables:

• •

$h_{i}(k)\in\mathbb{R}$: heat flow conveyed by ${\rm H}_{i}$ at time $k$ into room ${\rm Z}_{i}$,

• •

$\mathbf{h}(k)=[h_{1}(k)\ldots h_{m}(k)]^{\prime}\in\mathbb{R}^{m}$,

• •

$T_{i}(k)$: indoor temperature of zone ${\rm Z}_{i}$ at time $k$,

• •

$\mathbf{T}(k)=[T_{1}(k)\ldots T_{m}(k)]^{\prime}\in\mathbb{R}^{m}$,

• •

$[\underline{T}_{i}(k),\overline{T}_{i}(k)]$: thermal comfort range for ${\rm Z}_{i}$ at time $k$,

• •

$\underline{\mathbf{T}}(k)=[\underline{T}_{1}(k)\ldots\underline{T}_{m}(k)]^{\prime}\in\mathbb{R}^{m}$, $\overline{\mathbf{T}}(k)=[\overline{T}_{1}(k)\ldots\overline{T}_{m}(k)]^{\prime}\in\mathbb{R}^{m}$,

The dynamics of the indoor temperatures $\mathbf{T}(k)$ depend on the heat flows $\mathbf{h}(k)$ and on exogenous variables like outdoor temperature, solar radiation, appliances, lighting, occupancy, etc. For the sake of simplicity, available measurements/forecasts of some/all of such variables are collected in a vector $\mathbf{e}(k)$. Assuming linear time-invariant dynamics, the zone temperature vector evolution can be modeled in regressive form as:

[TABLE]

where the regression matrix $\mathbf{\Phi}(k)$ is given by

[TABLE]

being $k_{\mathbf{T}},k_{\mathbf{h}},k_{\mathbf{e}}$ suitable nonnegative integers that define the model order. The matrix $\Theta^{\mathbf{T}}$ collects the model parameters.

According to the above, thermal comfort achievement at time $k$ is expressed by the element-wise constraint

[TABLE]

It is appropriate to group constraint sets (1),(2),(3) as follows:

[TABLE]

3.2 HVAC system model

The HVAC system operates at constant water flow rate. With reference to Figure 1, let us define the following quantities:

• •

$T^{SND}(k)\in\mathbb{R}$: fluid temperature at the inlet of heat exchangers,

• •

$T^{RET}(k)\in\mathbb{R}$: temperature of return fluid,

• •

$v_{i}(k)\in\mathbb{R}$: heat exchanger actuation signal, i.e., commanded air flow at time $k$ for ${\rm H}_{i}$,

• •

$\mathbf{v}(k)=[v_{1}(k)\ldots v_{m}(k)]^{\prime}\in\mathbb{R}^{m}$,

• •

$\overline{v}_{i}$: maximum air flow rate allowed for ${\rm H}_{i}$, i.e.,

[TABLE]

For each zone Zi, the heat flow rate $h_{i}(k)$ can be expressed as

[TABLE]

where $\gamma_{i}$ is a coefficient of heating performance pertaining to heat exchanger Hi. Note that, in view of (5), enforcing limitation (4) for all zones yields the linear constraint sets

[TABLE]

and

[TABLE]

when the heat exchanger operates in heating and cooling mode, respectively.

3.2.1 Heating mode operation with thermal storage

The dynamics of the HVAC system in heating mode is now described. As pointed out in Section 2, the heat pump HPH is coupled to a thermal energy storage (TES). Define:

• •

$T^{HP_{H}}_{in}(k)\in\mathbb{R}$: thermal fluid temperature at HP inlet,

• •

$T^{HP_{H}}(k)\in\mathbb{R}$: thermal fluid temperature at HP outlet,

• •

$T^{HP_{H}}_{0}(k)\in\mathbb{R}$: HP outlet temperature reference, i.e., HP command signal,

• •

$W^{HP_{H}}(k)\in\mathbb{R}$: HP electrical energy consumption within the $k-$th time step,

• •

$T^{TES}(k)\in\mathbb{R}$: TES fluid temperature.

The inlet temperature of each heat exchanger Hi is assumed to be uniform and equal to the TES outlet temperature, i.e.,

[TABLE]

Under the realistic assumption that the time constants of the HP thermal fluid dynamics in response to a change in the reference $T^{HP}_{0}(k)$ are much smaller than the control sampling period $\tau_{s}$, the HP outlet fluid temperature can be expressed by

[TABLE]

so that the HP can be modeled as a unit time step delay system.

A further reasonable assumption is that the heat exchange at TES level is efficient enough such that

[TABLE]

It is worth noting that the above assumptions are indeed satisfied by the real-world devices emulated by EnergyPlus [26], the widely used realistic building simulation software used here for experimental validation.

The HP has to be switched on when the reference temperature is greater than the fluid temperature at the $HP$ inlet, i.e., as long as

[TABLE]

which leads to an electrical energy consumption

[TABLE]

where $\alpha^{HP_{H}}$ represents a device specification. Otherwise, the HP is switched off and $W^{HP_{H}}(k)=0$. Therefore, the condition

[TABLE]

can be assumed without loss of generality.

Concerning the TES dynamics, it is worth observing that $T^{TES}(k)$ depends on $T^{HP_{H}}(k)$ and on the return water temperature $T^{RET}(k)=T^{TES}_{in}$, which in turn dynamically depends on $T^{SND}(k)=T^{TES}(k)$ and on the total heat exchange at the heater level, i.e., $h(k)=h_{1}(k)+\dots+h_{m}(k)$. Therefore, the TES temperature dynamics can be modeled in regressive form as

[TABLE]

where the regression vector $\mathbf{\Phi}^{TES}(k)$ is given by

[TABLE]

being $\Theta^{TES}$ the (row) parameter vector and $k_{TT},~{}k_{TH},~{}k_{T\mathbf{h}}$ suitable model orders.

Given the models introduced above, and observing that the thermal energy stored in the TES is proportional to $T^{TES}(k)$, it is worth remarking that full control of individual zone temperatures and thermal energy storage can be achieved by manipulating the control variables $T^{HP_{H}}_{0}(k)$ and $\mathbf{v}(k)$.

3.2.2 Cooling mode operation

In order to derive a model of HVAC system in cooling mode, let us define the following:

• •

$T^{HP_{C}}_{in}(k)\in\mathbb{R}$: thermal fluid temperature at the inlet of the cooling heat pump,

• •

$T^{HP_{C}}(k)\in\mathbb{R}$: thermal fluid temperature at HP outlet,

• •

$T^{HP_{C}}_{0}(k)\in\mathbb{R}$: HP outlet temperature reference, i.e., HP command signal,

• •

$W^{HP_{C}}(k)\in\mathbb{R}$: HP electrical energy consumption within the $k-$th time step,

Similarly to the heating mode, the cooling HP dynamics can be modeled as:

[TABLE]

As long as $T^{HP_{C}}_{0}(k)\leq T^{HP_{C}}_{in}$, the HP is switched on and the corresponding electrical energy consumption is given by

[TABLE]

where $\alpha^{HP_{C}}$ pertains to the given device. Otherwise, it is switched off and $W^{HP_{C}}(k)=0$. Therefore, it can be assumed that

[TABLE]

In this mode, it holds that $T^{SND}(k)=T^{HP_{C}}(k)$ and $T^{RET}(k)=T^{HP_{C}}_{in}(k)$. The latter quantity dynamically depends on $T^{SND}(k)$ and on the total zone heat flow $h(k)=h_{1}(k)+\dots+h_{m}(k)$. Therefore, the following regressive model can be used for $T^{HP_{C}}_{in}(k)$:

[TABLE]

where the regression vector $\mathbf{\Phi}^{C}(k)$ is given by

[TABLE]

where $\Theta^{C}$ is the (row) parameter vector and $k_{TC},~{}k_{HC},~{}k_{C\mathbf{h}}$ define the model orders.

The constraint sets pertaining to the HVAC system operating in heating and cooling mode, respectively, can be grouped as follows:

[TABLE]

3.3 Electrical storage model

The characterization of the storage system based on rechargeable batteries used in this paper is as follows (see, e.g., [27, 28]). Let $E^{EES}(k)$ represent the state of charge at time $k$. Denote as $W^{EES}_{+}(k)\geq 0$ and $W^{EES}_{-}(k)\geq 0$ the battery charge and discharge signals, i.e., the amount of energy fed to and drawn from the battery in the $k-$th time interval, respectively. It holds that

[TABLE]

where $0<\eta<1$ is the battery efficiency. Furthermore, the amount of battery energy exchange in the same interval reads

[TABLE]

The following additional constraints

[TABLE]

[TABLE]

[TABLE]

where $\overline{W}^{ESS}_{+},\overline{W}^{ESS}_{-}$ and $\overline{E}^{EES}$ depend on the storage size and technology, and represent the maximum charge/discharge rate and capacity, respectively. The constraint sets pertaining to EES can be grouped as:

[TABLE]

In order for the model to be consistent, $W^{EES}_{+}(k)$ and $W^{EES}_{-}(k)$ cannot both be nonzero. This can be guaranteed by introducing the nonlinear constraint $W^{EES}_{+}(k)W^{EES}_{-}(k)=0$, $\forall k$. See Remark 2 in Section 4 on how the introduction of this constraint can be circumvented.

3.4 PV generation model

In this work, the well-known PVUSA model of a PV plant [29] is exploited. In this model, the maximum energy that can be drawn from the plant in the $k-$th time interval is expressed as a function of irradiance and environmental temperature as follows:

[TABLE]

where $I(k)$ is the global solar irradiance and $T^{A}(k)$ is the outside air temperature. Despite its simplicity, this model is very accurate when its parameters are fit to measured data and several efficient methods for their estimation under various conditions have been proposed (see [30],[31],[32] and references therein). Without loss of generality, it can be assumed that $I(k)$,$I^{2}(k)$, and $I(k)T_{a}(k)$ are part of the exogenous variable vector $\mathbf{e}(k)$. Therefore, (24) can be written in the form

[TABLE]

which represents a static linear-in-the-parameters model.

Denoting with $W^{PV}(k)$ the amount of energy drawn from the PV plant in the $k-$th interval, the following constraint must hold:

[TABLE]

The latter two equations are grouped in the constraint set

[TABLE]

3.5 Energy consumption and Demand-Response model

Let $W(k)$ represent the amount of energy drawn from the grid by the building in the $k-$th time interval. According to the model introduced above, the electrical energy balance equation reads:

[TABLE]

with $W(k)\geq 0$, since the possibility of injecting excess energy into the grid is not currently considered.

Let $p(k)$ represent the unit energy price or a forecast thereof, and consider a generic time horizon $\mathcal{I}(k,\lambda)$. The total energy drawn from the grid within $\mathcal{I}(k,\lambda)$ is given by

[TABLE]

and the total expected cost of energy in the same interval is equal to

[TABLE]

A Demand-Response model based on price-volume signals is considered in this paper. Such a model was introduced in [24] and is recalled next. A DR program is modeled as a sequence of DR requests ${\mathcal{R}}_{j}$, each carrying a time horizon $\mathcal{I}(h_{j},\mu_{j})$, a total energy bound $S_{j}$, and a monetary reward $R_{j}$. A request ${\mathcal{R}}_{j}$ is satisfied if the total building consumption within $\mathcal{I}(h_{j},\mu_{j})$, i.e., $W(h_{j},\mu_{j})$, is less or equal to the threshold $S_{j}$. In this case, a monetary reward $R_{j}$ is granted.

Definition 1

A DR program ${\mathcal{P}}$ is a sequence of DR requests ${\mathcal{R}}_{j}$, $j=1,2,\dots$, where ${\mathcal{R}}_{j}$ is the set

[TABLE]

being $\mathcal{I}(h_{j},\mu_{j})\subseteq\mathcal{K}$ and $\mathcal{I}(h_{j_{1}},\mu_{j_{1}})\cap\mathcal{I}(h_{j_{2}},\mu_{j_{2}})=\emptyset,~{}\forall j_{1}\neq j_{2}$.

The request ${\mathcal{R}}_{j}$ is satisfied if and only if

[TABLE]

For any given time horizon $\mathcal{I}(k,\lambda)$, let

[TABLE]

be the set of DR requests that occur within the time horizon. Moreover, define $\mathcal{J}(k,\lambda)$ as the set of indices identifying such DR requests, i.e.,

[TABLE]

To each request ${\mathcal{R}}_{j}$, let us associate a binary variable $\epsilon_{j}\in\mathcal{B}$ defined as

[TABLE]

The overall expected cost of operation of the building HVAC system within the time horizon $\mathcal{I}(k,\lambda)$ under the DR program $\mathcal{P}$, is therefore given by

[TABLE]

i.e., the expected cost of energy minus the total reward for the satisfied DR requests. Let $M$ be an a-priori known upper bound to the total building energy consumption per time step. Then, the set of constraints

[TABLE]

drives each variable $\epsilon_{j}$ to 1 if the respective DR request can be fulfilled, and to 0 otherwise.

The reader is referred to [24] for further details on this model and its MPC implementation.

4 Optimal HVAC and storage operation problem

The goal of this section is to devise an optimal schedule for the operation of the overall system (HVAC, PV, thermal and electrical storage) in order to minimize the building electricity bill under a DR program ${\mathcal{P}}$, while preserving comfort constraints. On a given interval $\mathcal{I}(k,\lambda)$, such a problem amounts to the optimal manipulation of the control variables $\mathbf{v}(l)$, $T^{HP_{H}}_{0}(l)$ $[T^{HP_{C}}_{0}(l)]$, $W^{EES}(l)$ for $l\in\mathcal{I}(k,\lambda)$ in order to minimize $C^{\mathcal{P}}(k,\lambda)$. To this purpose, it should be observed that the model is not linear in the above variables due to (5). Nevertheless, the model becomes linear if $\mathbf{h}(l)$ instead of $\mathbf{v}(l)$ are considered as decision variables. Once the optimal values $h_{i}^{*}(l)$ of $h_{i}(l)$ are available for zone Zi, the corresponding optimal heat exchanger actuation is given by

[TABLE]

In the formulation of the optimization problem, suitable forecasts or sensor measurements $\hat{\mathbf{e}}(k)$ of the exogenous variables are assumed to be available. Therefore, in the sequel, $\mathbf{e}(k)=\hat{\mathbf{e}}(k)$.

4.1 Heating mode operation with thermal storage

In order to formulate the optimal operation problem in heating mode over a time interval ${\cal I}(k,\lambda)$, the set of command variables $\mathbf{U}_{H}(k,\lambda)$, where

[TABLE]

is considered. Moreover, sensor measurements $\mathbf{T}_{m}(k)$, $T^{TES}_{m}(k)$, $T^{HP_{H}}_{m}(k)$, $E^{EES}_{m}(k)$ of $\mathbf{T}(k)$, $T^{TES}(k)$, $T^{HP_{H}}(k)$, $E^{EES}(k)$, respectively, are assumed to be available. Hence, the problem of minimizing the total energy cost $C^{\mathcal{P}}(k,\lambda)$ over the time interval $\mathcal{I}(k,\lambda)$ can be cast as a mixed-integer linear program (MILP) as follows.

Problem 1

Optimal combined heating/storage control under DR program ${\mathcal{P}}$.

[TABLE]

4.2 Cooling mode operation

In cooling mode, the set of control variables amounts to $\mathbf{U}_{C}(k,\lambda)$, where

[TABLE]

Sensor measurements $\mathbf{T}_{m}(k)$, $T^{HP_{C}}_{in,m}(k)$, $T^{HP_{C}}_{m}(k)$, $E^{EES}_{m}(k)$ of $\mathbf{T}(k)$, $T^{HP_{C}}_{in}(k)$, $T^{HP_{C}}(k)$, $E^{EES}(k)$, respectively, are considered. The problem of minimizing the total energy cost over $\mathcal{I}(k,\lambda)$ has the following formulation.

Problem 2

Optimal combined cooling/storage control under DR program ${\mathcal{P}}$.

[TABLE]

4.3 Managing constraint violations and uncertainty

The feasibility of the optimization problems introduced in the previous subsections cannot be guaranteed in general due to a number of circumstances that may occur in real use cases. Among such factors are the presence of uncertainties in the building component models, errors in sensor measurements or exogenous variable forecasts $\hat{\mathbf{e}}(k)$, inappropriate sizing or faults of the HVAC system, etc. Hence, a certain degree of constraint violation must always be accepted. In order to mitigate this issue, a common approach in MPC is to resort to soft constraints, that is, removing the relevant hard constraints and replace them with penalization terms in the objective function. However, in complex scenarios such as the one dealt with in this paper, designing suitable constraint penalization coefficients to be used in the cost function can be a very difficult task, and the optimal solution of the relaxed problem may turn out to be very sensitive to such design. To overcome this issue, a two-step procedure is adopted in this work, which is based on the observation that the constraints that may actually incur violations are the comfort constraints ${\cal C}(\cdot)$. In the first step, such constraints are relaxed by introducing a set of slack variables, and a modified version of Problem 1 or 2 is solved, in which the objective to be minimized is the $1-$norm of all constraint violations over ${\cal I}(k,\lambda)$. In the second step, Problem 1 or 2 is solved replacing the original comfort constraints with the optimal relaxed bounds computed in the first step. Note that this substitution always ensures feasibility of the latter problems.

In order to illustrate the procedure, the following set of positive slack variables pertaining to each time interval ${\cal I}(k,\lambda)$ is introduced:

[TABLE]

along with the relaxed constraint set

[TABLE]

Then, the following linear program, which is a modified version of Problems 1 and 2, is considered:

Problem 3

[TABLE]

where $\star=H$ for heating mode (Problem 1) and $\star=C$ for cooling mode (Problem 2). All other constraints (battery, DR, etc.) do not influence feasibility and can be omitted. Notice that the optimum of Problem 3 represents the minimum 1-norm of the constraint violations that must be accepted.

4.4 MPC algorithm

The proposed control procedure is implemented in a receding horizion fashion as the two-step optimization strategy in Algorithm 1.

Remark 1

At each step $k$, the proposed procedure consists of the solution of an LP (Problem 3) and of a MILP (Problem 1 or 2) involving a number of binary variables equal to the number of DR requests falling inside the prediction horizon ${\cal I}(k,\lambda)$. In realistic scenarios, the number of integer variables in the MILP is at most a few units [15] . The number of continuous variables as well as of constraints in all optimization problems scales in a linear fashion with the number of building zones. This makes the proposed procedure readily applicable to real-world scenarios involving large-scale buildings consisting of hundreds of zones, as demonstrated by the simulation experiments reported in Section 5.

Remark 2

The formulation of the optimization problems involved in Algorithm 1 does not take into account the constraint that the electrical storage cannot charge and discharge simultaneously, i.e., that $W^{EES}_{+}(l)$ and $W^{EES}_{-}(l)$ cannot be both nonzero. However, unless $E^{EES}(l+1)=0$ or $E^{EES}(l+1)=\overline{E}^{EES}$, such condition is always met because the battery efficiency satisfies $\eta<1$. On the contrary, when the battery state of charge hits one of the bounds, $W^{EES}(l)$ can be set such that $E^{EES}(l+1)=0$ or $E^{EES}(l+1)=\overline{E}^{EES}$ according to (19),(20).

5 Test cases

In this section, the application of the proposed control technique to a large-scale building is presented. Heating and cooling operations are considered in Sections 5.1 and 5.2, respectively, while Section 5.3 is devoted to analyzing the performance of the proposed method when inaccurate forecasts of exogenous variable are available. A discussion of the results is given in Section 6.

The simulated test structure is an 8-floor building divided in 126 zones (Fig. 2) located in Turin, Italy. The first two floors are devoted to commercial and office activities, respectively, while the remaining floors are assigned to residential apartments. Plans of the three types of floors are reported in Fig. 3.

The building is equipped with a PV plant, an electrical storage system and two heat pumps used for heating and cooling purposes, respectively. The heating system is also equipped with a TES to enable hot water storage. Heating and cooling circuit schemes are shown in Fig. 4 and 5. The building has been designed using DesignBuilder [33] while simulations have been performed via EnergyPlus [26] connected to Matlab through the Building Controls Virtual Test Bed (BCVTB) [34]. Building features are reported in Table 1, while materials are summarized in Table 2. The energy price time series $p(k)$ has been taken from the Italian electricity market [35]. The EnergyPlus simulation model is assumed as the real building. The sampling time $\tau_{s}$ is set to 10 minutes, while the horizon length used by the MPC is fixed to 12 hours, corresponding to 72 samples, i.e., $\lambda=72$ in Algorithm 1.

5.1 Winter season simulation - heating operation

The thermal behavior of the building has been modeled as in Section 3.1 by means of the autoregressive model (1)-(2). To this purpose, it is appropriate to define the exogenous input vector $\mathbf{e}(k)$ as

[TABLE]

where $T^{A}(k)$ is the outside air temperature, $I(k)$ is the solar irradiance and $\mathbf{G}(k)$ denotes the vector of internal heat gains of the building zones due to human occupancy, lighting and equipments.

For this model, an identification experiment has been performed over a time horizon of 18 days, which have been split in 14 days for estimation and 4 for validation. For each building zone, a coupled model involving all neighboring zones has been employed. The standard Best FIT index (FIT) defined in [36, 37] has been used to assess the quality of the estimated model. The average value of the FIT index for all the identified building zones is 75%, 68% and 66%, for 1, 6 and 12-hour ahead predictions, respectively. For a qualitative evaluation, in Fig. 6, the 6-hour ahead prediction of the internal temperature of an office zone is compared with the real behavior during the model validation phase.

Concerning the TES, model (12)-(13) has been identified, and the associated FIT index turns out to be over 90% even for 12-hour ahead prediction. The PV plant has been identified by employing model (25).

Depending on the type of zone (commercial, office, residential), different comfort constraints are enforced. For commercial and office zones, temperature bounds are set to 20-24°C from 8:00 to 18:00, and 15-24°C elsewhen. Bounds for residential zones are 20-24°C from 7:00 to 9:00 and from 19:00 to 01:00, and 15-24°C otherwise. It is worthwhile to note that such bounds can be freely adjusted for each zone to adapt to real scenarios.

A 3-day simulation has been performed to evaluate the proposed control strategy. The price/volume DR requests in the considered period are assumed to be known one day in advance [15] and are reported in Table 3. The identified model has been used by the MPC to generate fan input signals, HP setpoints and EES charging/discharging commands, while real measured data have been generated by the EnergyPlus simulation of the designed building.

At this stage, perfect knowledge of exogenous inputs is assumed. Such an assumption, although not realistic, has been enforced to better assess the quality of the proposed technique. Simulations in presence of uncertain forecasts will be analyzed in detail in Section 5.3.

In Fig. 7, outdoor temperature $T^{A}$ and solar irradiance $I$ for the considered days are reported, while in Fig. 8, TES temperature $T^{TES}$, HP setpoint $T^{HP_{H}}_{0}$, EES state of charge $E^{EES}$, energy drawn from the PV plant $W^{PV}$, unit energy price $p$ and total energy drawn from the grid $W$ are shown. In Fig. 9, the temperature profiles $T_{i}(k)$ of three different sample zones are depicted along with the required comfort bounds $[\underline{T}_{i},\overline{T}_{i}]$ and fan actuation $v_{i}$. It can be noted that the zone temperatures lie within the comfort bounds with small deviations. Indeed, the average bound violation of the worst performing zone over three days is about 0.13°C, see Table 4. Notice that, since exact input forecasts are assumed, such infringements are due to modeling errors. This demonstrates the ability of the bound relaxation technique in Section 4.3 to manage unpredictable bound violations by pushing the zone temperatures back within the prescribed bounds in a short time after the violations occur. It is worthwhile to recall that keeping air temperature within hard comfort bounds at all times is unfeasible for any conceivable control strategy due to modeling errors and/or unpredictable disturbances. In Table 4, numerical values of total cost, DR rewards, and comfort bound violations are reported.

To further assess the performance of the proposed procedure, a standard thermostatic controller with 0.5°C hysteresis has been used as a benchmark. Since thermostatic control is clearly unable to track step variations of the comfort constraints, preheating/precooling is performed with suitable advance in order to provide a fair comparison. Moreover, when thermostatic control is in place, the EES is operated in the following fashion: if the power currently produced by the PV panels exceeds that needed by the HP, the surplus is stored in the battery, otherwise the controller draws power from the battery in the first place, and then from the grid. Since DR requests cannot be handled by thermostatic control, simulations in absence of DR programs have been performed for the sake of fairness.

In Table 5, the total cost is summarized for both strategies along with the average bound violation of the worst performing zone. In the considered days, the total cost achieved by adopting the proposed control technique turns out to be 25.49% less than that obtained by the benchmark.

5.2 Summer season simulation - cooling operation

As described in Section 3.2.2, the problem structure for cooling mode is similar to the heating case, with the exception that no TES is assumed to be present. The identified models used under this condition have been obtained in a similar manner. A model for the dynamics of the fluid temperature at the HP inlet ($T^{HP_{C}}_{in}$) is identified according to (17)-(18), and the resulting FIT index turns out to be over 80% for 12-hour ahead predictions. This step is not explicitly needed for the heating case since, according to (8), the fluid temperature at HP inlet is the same as the TES temperature $T^{TES}$.

For commercial and office zones, comfort bounds are set to 22-24°C from 8:00 to 18:00, and 22-28°C elsewhen. Bounds for residential zones are 22-24°C from 7:00 to 9:00 and from 19:00 to 01:00, and 22-28°C otherwise.

In Table 6 and Fig. 10-11, the results of a three-day simulation are summarized. In Table 7, the comparison with the thermostatic controller is shown. As for the winter season, the proposed MPC reduces the cost by 35.67% with respect to the benchmark, still maintaining similar zone comfort.

5.3 Performance analysis under uncertainties

Up to this point, exact forecasts of the exogenous inputs $\mathbf{e}$ have been considered. Such an assumption is not acceptable in a realistic scenario for weather variables and zone internal gains. On the contrary, it is reasonable to assume that energy price is exactly known one day in advance, as well as the daily DR program [15]. For these reasons, simulations have been performed to evaluate how the proposed technique is sensitive to forecasting errors in outdoor temperature, solar irradiance and internal gains.

In a real scenario, predictions of external temperature are provided by national weather forecast services. Outdoor temperature forecasting errors have been simulated as follows. First, a second-order autoregressive model has been used to compute the noise signal

[TABLE]

where $\alpha_{1}$ and $\alpha_{2}$ are the model coefficients and ${\varepsilon}$ is a zero-mean Gaussian distributed random variable. Notice that it is quite common to use autoregressive systems to model temperature uncertainty, see e.g., [38]. Then, denoting by $T_{true}^{A}$ the real temperature profile, the temperature forecast $\widehat{T}^{A}$ to be used by the MPC has been obtained as

[TABLE]

Notice that $\widehat{T}^{A}(k)=T_{true}^{A}(k)$, since the outdoor temperature is measured. On the other hand, the term $\frac{l-k}{\lambda}$ takes into account the fact that forecasting accuracy decreases with time. State-of-the-art models for temperature forecasting are able to generate day-ahead estimates with an error of about 2°C [39]. Hence, the coefficients in (60) have been chosen conservatively in order obtain a maximum error of 3°C over a 12-hour ahead prediction.

To test the behavior of the proposed control law, a number of uncertain input profiles have been generated and used in simulations. In Fig. 12, five different 12-hour ahead forecasts are reported, along with the true temperature profile.

Predictions of solar irradiance have been obtained by multiplying the real irradiance by a signal generated as in (60). In Fig. 13, different profiles of solar irradiance predictions are reported along with the actual signal.

To obtain forecasts on zone internal gains, a noise signal as in (60) has been added to the real internal gain, by saturating the perturbed signal to 0. In Fig. 14, the true internal gain is depicted for each zone type along with five forecasts.

Twenty simulations in heating mode have been generated for each uncertain input. In Table 8, the worst-case value of the total cost and of the average comfort bound violation of the worst performing zone are reported. Inaccurate forecasts are considered both individually and combined. The box plots related to the overall cost for the considered scenarios are shown in Fig. 15.

It can be noticed that the presence of inaccurate forecasts may raise the overall cost up to 7.79% with respect to the nominal case, while the difference on comfort is negligible. The latter fact shows the ability of the control algorithm to minimize the thermal discomfort when unexpected situations occur.

6 Discussion

In this section, the results of the simulation experiments are analyzed and discussed.

Let us first consider the simulation reported in Section 5.1, corresponding to winter season and perfect knowledge of exogenous inputs. The identified models are able to predict the system dynamics accurately. As expected, the performance indexes decrease as the prediction horizon grows, yet they remain on acceptable values. Moreover, it is worth remembering that in a receding horizon approach, predictions are updated at each time step, thus making it possible for the control algorithm to adapt against inaccuracies on long-term forecasts.

Concerning the control system simulation in Fig. 8, it can be observed that the EES and the TES are charged and discharged twice a day in order to take advantage of electricity price fluctuations. Moreover, it is apparent that both the EES and the TES play an important role in order to satisfy DR requests. From Fig. 8 and 9, it is apparent that the last DR request is not fulfilled, meaning that consumption reduction during the DR interval is not deemed profitable by the controller. In Fig. 9, one can see that the comfort bounds are mostly satisfied. In this case, since exact forecasts are assumed, bound violations are only due to discrepancies between the real system and the identified model. However, the magnitude of such bound violations is almost negligible. In fact, as reported in Table 4, the average bound violation for the worst zone is less than 0.14°C. The proposed control technique leads to a total cost which is 25.49% less than that obtained by standard thermostatic rules, as reported in Table 5. Since the benchmark cannot exploit the knowledge of DR requests, no DR program has been assumed in this case to ensure a fair comparison. Of course, the gap between the two strategies in the presence of DR programs is larger.

Since the assumption of exact forecasts is unrealistic, the cost associated to the benchmark controller is also compared with that obtained by the MPC under inaccurate input predictions. To this purpose, Table 8 and Fig. 15 show that the MPC cost under uncertain forecasts does not change significantly w.r.t. the nominal one. In fact, the maximum cost under inaccurate forecasts is less than 8% greater than that obtained for exact forecasts. This is a consequence of the receding horizon strategy, which implies that only the prediction errors occurring in the near future significantly affect the performance. This makes the proposed approach intrinsically robust to uncertainty, since short-term forecasts are usually quite accurate.

In addition to affecting the total cost, uncertainty also plays a role in the preservation of zone temperature comfort bounds. However, this effect is mitigated by the MPC implementation described in Section 4.3, as it is clear from Table 8, which shows that only a slight increment of discomfort occurs in the presence of unreliable forecasts.

Similar considerations can be made for the results obtained in cooling operation, i.e., in summer season.

A final remark concerns the computational burden associated with the controller implementation. For the considered 126-zone building, the time required for computations at each step (i.e., every 10 minutes) is about 20 seconds on a standard PC111Computations have been performed using CPLEX [40] to solve the LPs, on an Intel Core i7-3770 at 3.40 GHz with 8 GB of RAM., showing that the proposed algorithm can be efficiently adopted in large-scale applications.

7 Conclusion and future research

In this paper, the problem of optimizing the operation cost (i.e., the electrical energy bill) of a building integrating a centralized HVAC, thermal and electrical storage facilities, and PV generation has been addressed. Participation in a DR program has also been considered. The proposed approach exploits a receding horizon control MPC strategy involving at each step the solution of an LP and of a MILP with a number of integer variables equal to the number of DR requests in the prediction horizon. The procedure can be fruitfully applied to large size buildings. Experimental validation using a realistic simulation framework has been carried out.

Further studies will address more complex setups in which the building is considered as a microgrid in its own right, including electric vehicles, appliances and other kinds of loads. Special attention to active/reactive power flow should be paid in this case in order to guarantee the satisfaction of electrical constraints. More general DR scenarios will also be addressed, such as the presence of incentives for keeping energy consumption above a minimum, which can be efficiently achieved using storage facilities. Finally, different methods for handling uncertainty, like chance constraints or scenario-based approaches, will be investigated and compared with the method proposed in this paper. To this purpose, suitable problem formulations and related relaxations will be studied to enable such techniques to manage large-scale buildings.