## Abstract

This paper presents a novel architecture for model predictive control (MPC)-based indoor climate control of multi-zone buildings to provide energy efficiency. Unlike prior works, we do not assume the availability of a high-resolution multi-zone building model, which is challenging to obtain. Instead, the architecture uses a low-resolution model of the building that is divided into a small number of “meta-zones” that can be easily identified using existing data-driven modeling techniques. The proposed architecture is hierarchical. At the higher level, an MPC controller uses the low-resolution model to make decisions for the air handling unit (AHU) and the meta-zones. Since the meta-zones are fictitious, a lower level controller converts the high-level MPC decisions into commands for the individual zones by solving a projection problem that strikes a trade-off between two potentially conflicting goals: the AHU-level decisions made by the MPC are respected while the climate of the individual zones is maintained within the comfort bounds. The performance of the proposed controller is assessed via simulations in a high-fidelity simulation testbed and compared to that of a rule-based controller that is used in practice. Simulations in multiple weather conditions show the effectiveness of the proposed controller in terms of energy savings, climate control, and computational tractability.

## 1 Introduction

The application of model predictive control (MPC) for commercial heating, ventilation, and air conditioning (HVAC) systems for both energy efficiency and demand flexibility has been an active area of research; see the review articles [1,2] and the references therein. Our paper contributes to this literature, specifically to that on model predictive control of air handling units (AHUs) and the zones they serve through variable air volume (VAV) boxes. Figure 1 shows an illustration of such a system, which are common in large commercial buildings.

Fig. 1
Fig. 1
Close modal

Several of the MPC formulations proposed in the past are for buildings with one zone [36] or a small number of zones [7,8]. A direct extension of such formulations for large multi-zone buildings has two main challenges. First, solving the underlying optimization problem in MPC for a building with a large number of zones is computationally complex because of the large number of decision variables. To reduce the computational complexity, several distributed and hierarchical approaches have been proposed [915]. The second challenge, which has attracted far less attention, is that MPC requires a “high-resolution” model of the thermal dynamics of a multi-zone building. High-resolution means that the temperature of every zone in the building is a state in the model and the control commands for every zone are inputs in the model. One way of obtaining such a multi-zone model is by first constructing a “white box” model, such as by using a building energy modeling software, and then simplifying it to make it suitable for MPC, e.g., Ref. [16]. But constructing a white box model is expensive; it requires significant effort [17]. Moreover, the resulting model may not reflect the building as is. Another way of obtaining a high-resolution multi-zone model is by utilizing data-driven techniques, which use input-output measurements. Getting reliable estimates using data-driven modeling is challenging even for a single-zone building, as a building’s thermal dynamics is affected by a non-trivial and unmeasurable disturbance, the heat gain from occupants and their use of equipment, that strongly affects quality of the identified model [18]. In the case of multi-zone model identification, it becomes intractable since the model has too many degrees-of-freedom: as many unknown disturbance signals as there are number of zones. To the best of our knowledge, there are no works on reliable identification of multi-zone building models without making assumptions on the nature of the disturbance affecting individual zones [19].

In addition to the challenges mentioned above, most of the prior works—whether on single-zone or on multi-zone buildings—ignore humidity and latent heat in their MPC formulations. The inclusion of moisture requires a computationally convenient cooling and dehumidifying coil model. MPC formulations that exclude humidity can lead to poor humidity control, or higher energy usage as they are unaware of the latent load on the cooling coil [5].

In this work, we propose a humidity- and latent heat-aware MPC formulation for a multi-zone building with a VAV HVAC system, such as the one in Fig. 1. Cooling and dehumidification of air occur only at the AHU, using chilled water. Heating of air occurs only at the VAV boxes and is provided through either heating hot water or electricity. The problem we consider is making optimal decisions in real-time for actuators in both the AHU and the VAV boxes; details of control commands are described in Sec. 2.

A two-level control architecture is proposed here, consisting of a high-level controller (HLC) and a low-level controller (LLC), to address the challenges mentioned earlier. The high-level controller decides the control commands for air handling unit, while the low-level controller decides the control commands for individual zones. The HLC is an MPC controller that uses a “low-resolution” model of the building with a small number of “meta-zones”, with each meta-zone being a single-zone equivalent of a part of the building consisting of several zones. For instance, in the case study presented here, a 33-zone three-floor building is aggregated to a three meta-zone model, with each meta-zone corresponding to a floor. The advantage of such an approach is that a high-resolution multi-zone model is not needed as a starting point. Rather, a single-zone equivalent model of each meta-zone, in which disturbance in all the zones are aggregated into one signal, can be directly identified from measurements collected from the building. The identification problem of such a single-zone equivalent model is more tractable [20]. In this paper, we use the system identification method from Ref. [20], though other identification methods can also be used. Since the HLC uses a low-resolution model with a much smaller number of meta-zones than that in the building, its computational complexity is low. However, this reduction of computational complexity creates a different challenge. Since the decision variables of the optimization problem in the HLC correspond to the meta-zones (air flowrate, temperature, etc.), they do not correspond to those for the actual zones of the building. The LLC resolves this difficulty by solving a projection problem that appropriately distributes aggregate quantities computed by the HLC to commands for individual zones. The LLC uses feedback from each zone to assess their needs and ensures that indoor climate of each zone is maintained and the aggregate quantities for each meta-zone are consistent with what the HLC computed for that meta-zone.

The proposed controller—that includes the HLC and LLC—is hereafter referred to as MZHC which stands for multi-zone hierarchical controller. Its performance is assessed through simulations on a “virtual building” plant. The plant is representative of a section of the Innovation Hub building comprising of 33 zones and is located at the University of Florida campus. The plant is constructed using modelica [21]. The performance of the proposed controller is compared with that of Dual Maximum controller as a baseline [22]. The dual maximum controller—which is referred to as BL (for baseline)—is a rule-based controller, and is one of the more energy efficient controllers among those used in practice [22]. Simulation results show that using the proposed controller provides significant energy savings when compared to BL while maintaining indoor climate.

Compared to the literature on MPC design for multi-zone building HVAC systems, our work makes three principal contributions, with details discussed in Sec. 1.1. (i) The first contribution is that the proposed method does not assume availability of a high-resolution model of the multi-zone building which is difficult to obtain. Instead, it can utilize existing data-driven methods that can quickly identify a low-resolution model of the multi-zone building from measurements. (ii) Since the MPC part of the proposed controller uses a low-resolution model with a small number of meta-zones, the method is scalable to buildings with a large number of zones. Although distributed iterative computation has been proposed in the literature as an alternate approach to reducing computational complexity, ours can be solved in a centralized setting. (iii) The third contribution is the incorporation of humidity and latent heat in our multi-zone MPC formulation, which has been largely ignored in the literature on MPC for buildings, and especially so in the literature on multi-zone building MPC. Our simulations show that when using MZHC, the indoor humidity constraint is active, especially during hot-humid weather. Without humidity being explicitly considered, the controller is likely to have caused high space humidity in an effort to reduce energy use.

The rest of this paper is organized as follows. Section 1.1 discusses our work in relation to the literature on multi-zone MPC. Section 2 describes a multi-zone building equipped with a VAV HVAC system and the models we use in simulating the plant (the system to be controlled). Section 3 presents the proposed MPC-based hierarchical controller. Section 4 describes a rule-based baseline controller with which the performance of the proposed controller is compared. The simulation setup is described in Sec. 5. Simulation results are presented and discussed in Sec. 6. Finally, the main conclusions are provided in Sec. 7.

### 1.1 Comparison With Literature on Multi-Zone Model Predictive Control.

The works on MPC for multi-zone buildings can be roughly categorized into three bins based on the control commands considered in their formulation: zone-level commands only, AHU-level commands only, and both AHU- and zone-level commands.

The first category of works use MPC to vary only the zone-level control commands such as the zone temperature set points and supply airflow rates but not the AHU-level commands such as conditioned air temperature and outside airflow rate [7,9,10,1215,23]. Since these works do not consider all the available degrees-of-freedom, their energy savings potential is lower.

The second category is a complement of the first: MPC is used to vary only the AHU-level control inputs but not the zone-level control inputs. The controller presented in Ref. [24] belongs to this category. Focusing only on AHU-level inputs helps them avoid the need to identify a high-resolution multi-zone model. However, excluding zone-level control inputs reduce energy savings potential since the controller cannot use all available degrees-of-freedom.

The third category of papers, which this work belongs to, considers both zone-level and AHU-level control inputs in their formulation [8,11,25]. In Refs. [8,25], centralized approaches are used as multi-zone buildings with only a small number of zones are considered. In Ref. [11], a distributed algorithm is proposed to reduce the computational complexity of multi-zone building MPC, but their algorithm requires a high-resolution multi-zone model.

Besides, none of the papers mentioned earlier take into account humidity and latent heat explicitly in their optimization problem formulation.

## 2 System and Problem Description, and Plant Simulator

Our focus is a multi-zone building equipped with a VAV HVAC system, whose schematic is shown in Fig. 1. In such a system, part of the air exhausted from the zones is recirculated and mixed with outdoor air. This mixed air is sent through the cooling coil where the air is cooled and dehumidified to the conditioned air temperature (Tca) and humidity ratio (Wca). This conditioned air is then sent through the supply ducts to the VAV boxes, which have a damper to control airflow, and finally supplied to the zones. Some VAV boxes have reheat coils; they can change temperature of supply air but not humidity, i.e., Tsa,iTca and Wsa,i = Wca, where Tsa,i and Wsa,i are the temperature and humidity ratio of supply air to the ith zone. If a VAV box is not equipped with a reheat coil (cooling only), then the temperature of air supplied by it to its zone will be at the conditioned air temperature, i.e., Tsa,i = Tca.

The control commands for a multi-zone VAV HVAC system with nz zones (i.e., VAV boxes) are as follows:
$u:=(moa,Tca,msa,i,Tsa,i,i=1,…,nz)$
(1)
where moa is the outdoor airflow rate, Tca is the conditioned air temperature, msa,i is the supply airflow rate to the ith zone, and Tsa,i is the supply air temperature to the ith zone. Note that the humidity of conditioned air (Wca) which is supplied to all the zones is indirectly controlled through Tca. Of the nz VAVs/zones in the building, $nzrh$ VAVs are equipped with a reheat coil and $nz−nzrh$ VAVs do not have a reheat coil (cooling only). For the latter, the supply air temperature will be the same as the conditioned air temperature, i.e., Tsa,i = Tca.

The control commands in Eq. (1) are sent as set points to the low-level control loops which are typically comprised of proportional integral (PI) controllers. Control of the supply airflow rates msa,i, i = 1, …, nz, is performed through coordination between VAV boxes and the fan-speed PI controller at the AHU. When the dampers of a VAV box are opened by its local controller, the duct static pressure pduct decreases (see Fig. 1), and the fan-speed controller increases the fan’s speed to maintain the duct static pressure at its setpoint, in effect increasing the air flowrate. The PI controllers at the VAV boxes thus track their local supply airflow rate setpoints. The outdoor airflow rate moa is controlled by another PI loop that varies the position of outdoor air and return air dampers to maintain this setpoint. The conditioned air temperature Tca is controlled by a PI loop that varies the chilled water valve position to maintain this setpoint. The supply air temperatures $Tsa,i,$$i=1,…,nzrh$ are controlled by local PI loops at the VAV boxes that vary the reheat valve positions (or reheat powers).

The role of a climate control system is to vary the control commands in Eq. (1) so that three important and somewhat conflicting goals are satisfied: (i) ensure thermal comfort, (ii) maintain indoor air quality, and (iii) use minimum amount of energy/cost.

The paper uses a discrete time MPC setting, with k denoting discrete time. The time interval between two consecutive discrete time indices k and k + 1 is denoted by Δt and is called the model interval. The control commands are updated less frequently, and the time interval between two discrete time indices in which the control commands are updated is denoted by ΔT, an integer multiple of Δt, and is called the control interval [26].

### 2.1 Virtual Building (Simulator).

The virtual building (VB) is a high-fidelity model of a building’s thermal dynamics and its HVAC system that will act as the plant for the controllers. The VB is chosen to mimic part of the Innovation Hub building in Gainesville, FL, which is serviced by AHU-2 (among the two AHUs that serve Phase I). Figures 2 and 3 show photos of the building and the relevant floor plans, respectively. The rooms supplied by the same VAV box are grouped together to form one large space (zone); the zones are enclosed by dotted lines in Fig. 3. The first floor has 15 rooms that are grouped into nine zones, the second floor has 20 rooms which are grouped into 12 zones, and the third floor has 21 rooms which are grouped into 12 zones. In total, there are 56 rooms grouped into 33 zones. The virtual building thus consists of an air handling unit and 33 VAV boxes, of which 29 are equipped with reheat coils, and the remaining four do not have reheat coils (cooling only). The zones primarily consist of offices and labs.

Fig. 2
Fig. 2
Close modal
Fig. 3
Fig. 3
Close modal

We use the modelica library IDEAS (Integrated District Energy Assessment by Simulation) [27] to model the building’s thermal dynamics. The core of IDEAS library is the modelica IBPSA library [28] which is also used by other modelica building libraries such as the Buildings library [29] from LBNL (Lawrence Berkeley National Laboratory). We used the IDEAS library in this work since we found its graphical user interface (GUI) to be intuitive and simple to use to create a model which mimics the configuration of the Innovation Hub building.

Inputs to the building thermal dynamics portion of the VB are supply airflow rate (msa,i), supply air temperature (Tsa,i), and supply air humidity ratio (Wsa,i) for all the zones, and the exogenous disturbances, which are (i) solar irradiation (ηsol), outdoor air temperature (Toa), and outdoor air humidity ratio (Woa), (ii) internal sensible and latent heat loads due to occupants, which are computed based on the number of occupants provided to the zone block, and (iii) internal heat load due to lighting and equipment. Outputs of the simulator are temperature (Tz,i) and humidity ratio (Wz,i) of all the zones. Figure 4 shows floor 1 of the VB created in dymola using components from the IDEAS library [27]. We omit the details of the VB construction using modelica/dymola; the interested reader is referred to Ref. [30].

Fig. 4
Fig. 4
Close modal

#### a Cooling and Dehumidifying Coil Model.

The cooling coil model has five inputs and two outputs. The inputs are supply airflow rate (msa), mixed air temperature (Tma), humidity ratio (Wma), chilled water flowrate (mw), and chilled water inlet temperature (Twi). The outputs are conditioned air temperature (Tca) and humidity ratio (Wca). We use a gray box data-driven model developed in our prior work [5]. The interested readers are referred to Sec. 2.1.2 of Ref. [5] for details regarding the model.

#### b Power Consumption Models.

For the HVAC system configuration presented in Fig. 1, there are three main components which consume power. They are supply fan, cooling and dehumidifying coil, and reheating coils. The fan power consumption is modeled as follows:
$Pfan(k)=αfanmsa(k)3$
(2)
where msa(k) is the total supply airflow rate at the AHU [31].
The cooling and dehumidifying coil power consumption is modeled to be proportional to the heat it extracts from the mixed air stream:
$Pcc(k)=msa(k)[hma(k)−hca(k)]ηccCOPc$
(3)
where hma(k) and hca(k) are the specific enthalpies of the mixed and conditioned air, respectively; ηcc is the cooling coil efficiency, and COPc is the chiller coefficient of performance. Since a part of the return air is mixed with the outside air, the specific enthalpy of the mixed air is as follows:
$hma(k)=roa(k)hoa(k)+(1−roa(k))hra(k)$
(4)
where hoa(k) and hra(k) are the specific enthalpies of the outdoor and return air, respectively, and roa(k) is the outside air ratio: $roa(k):=moa(k)/msa(k)$. The specific enthalpy of moist air with temperature T and humidity ratio W is given by Ref. [32]: $h(T,W)=CpaT+W(gH2O+CpwT)$, where $gH2O$ is the heat of evaporation of water at 0 °C, and Cpa, Cpw are specific heat of air and water at constant pressure.
The reheating coil power consumption in the ith VAV box is modeled to be proportional to the heat it adds to the conditioned air stream:
$Preheat,i(k)=msa,i(k)Cpa[Tsa,i(k)−Tca(k)]ηreheatCOPh$
(5)
where ηreheat is the reheating coil efficiency and COPh is the boiler coefficient of performance.

#### c Overall Plant.

The overall plant, i.e., virtual building—consisting of the building thermal model, cooling and dehumidifying coil model, and power consumption models—is simulated using simulink and matlab. The building thermal model is constructed in dymola 2021 and it is exported into an FMU (Functional Mockup Unit). It is then imported into simulink using the FMI Kit for simulink. The remaining models are constructed directly in simulink.

## 3 Proposed Multi-Zone Hierarchical Control

Recall that both the proposed and the baseline controllers need to decide the following control commands:
$u(k):=[moa(k),Tca(k),msa,i(k),Tsa,i(k)]T∈ℜ2+nz+nzrh$
Figure 5 shows the structure of the proposed MZHC. The high-level controller is based on MPC and decides the control commands for the AHU: outdoor air flowrate (moa) and conditioned air temperature (Tca). The low-level controller is a projection-based feedback controller and decides the control commands for each of the VAV boxes/zones: supply air flowrate (msa,i) and supply air temperature (Tsa,i). These controllers are described in detail next.
Fig. 5
Fig. 5
Close modal

### 3.1 Model Predictive Control-Based High-Level Controller.

The HLC is based on MPC that uses a low-resolution model of the building which is divided into a small number of meta-zones. Each meta-zone is an aggregation of multiple zones in the real building. This aggregation can be done in any number of ways. In this work, we aggregate all the zones in a floor into a meta-zone, which is denoted by $f∈F:={1,…,nf}$, where nf is the total number of floors/meta-zones. The set of all VAVs/zones in floor f is denoted as If (so $|∪f∈FIf|=nz$), of which those equipped with reheat coils is denoted as Irh,f (so $|∪f∈FIrh,f|=nzrh$). The HLC decides on the following control commands based on the aggregate models:
$uHLC(k):=(moa(k),Tca(k),msa,f(k),Tsa,f(k),∀f∈F)∈ℜ2+(2×nf)$
(6)
where $msa,f(k):=∑i∈Ifmsa,i(k)$ is the aggregate (total) supply airflow rate to all the zones in floor/meta-zone f and Tsa,f(k) is the aggregate supply air temperature. Of the control commands computed in Eq. (6), moa(k) and Tca(k) can be directly sent to the plant. The remaining information computed by the HLC including msa,f(k) and Tsa,f(k) are used by the LLC, described in Sec. 3.2, to decide on the supply airflow rate (msa,i(k)) and supply air temperature (Tsa,i(k)) for the individual zones/VAV boxes in each floor.

A comment on notation: all variables with a subscript i are for the individual zones, while the variables with a subscript f represent the aggregate quantities for each meta-zone.

For the MPC formulation, we use a model interval of Δt = 5 min, a control interval of ΔT = 15 min, and a prediction/planning interval of T = 24 h. So, we have T = NΔT and ΔT = MΔt, where N = 96 (planning horizon) and M = 3. The control inputs for N time-steps are obtained by solving an optimization problem of minimizing the energy consumption subject to thermal comfort, indoor air quality, and actuator constraints. Then, the control commands obtained for the first time-step are sent to the plant and the LLC. The optimization problem is solved again for the next N time-steps with the initial states of the model obtained from a state estimator, which uses measurements from the plant. This process is repeated at the next control time-step, i.e., after an interval of ΔT.

To describe the optimization problem, first we define the state vector x(k) and the vector of control commands and internal variables v(k) as follows:
$x(k):=(Tz,f(k),Tw,f(k),Wz,f(k),∀f∈F)∈ℜ3×nf$
(7)
$v(k):=(uHLC(k),mw(k),Wca(k))∈ℜ2+(2×nf)+2$
(8)
where Tz,f(k), Tw,f(k), and Wz,f(k) are the aggregate zone temperature, wall temperature, and humidity ratio of floor/meta-zone f, respectively; uHLC(k) is the control command vector defined in Eq. (6) and mw(k) is the chilled water flowrate into the cooling coil. The exogenous input vector is as follows:
$w(k):=(ηsol(k),Toa(k),Woa(k),qint,f(k),ωint,f(k),∀f∈F)∈ℜ3+(2×nf)$
(9)
where ηsol(k) is the solar irradiance, Toa(k) is the outdoor air temperature, Woa(k) is the outdoor air humidity ratio, qint,f(k) is the aggregate internal heat load in floor/meta-zone f due to occupants, lights, equipment, etc., and ωint,f(k) is the aggregate rate of water vapor generation in floor/meta-zone f due to occupants and other sources. We denote the forecast of these exogenous inputs as $w^^$; in Sec. 5, we discuss how these forecasts are obtained. The vector of nonnegative slack variables $ζ(k):=(ζT,flow(k),ζT,fhigh(k),ζW,flow(k),$$ζW,fhigh(k),∀f∈F)∈ℜ4×nf,$ is introduced for feasibility of the optimization problem.
The optimization problem at time index j is the following minimization problem subject to a number of constraints that will be described later:
$minV,X,Z∑k=jj+NM−1[Pfan(k)+Pcc(k)+∑f∈FPreheat,f(k)]Δt+Pslack(k)$
(10a)
where Pfan(k) is given by Eq. (2), Pcc(k) is given by Eq. (3), $Preheat,f(k):=msa,f(k)Cpa[Tsa,f(k)−Tca(k)]/ηreheatCOPh,V:=[vT(j),$$vT(j+1),…,vT(j+NM−1)]T,X:=[xT(j+1),$$xT(j+2),…,xT(j+NM)]T,$ and $Z:=[ζT(j+1),ζT(j+2),…ζT(j+NM)]T$. The last term, Pslack, penalizes the aggregate zone temperature and humidity slack variables:
$Pslack(k):=∑f∈F[λTlowζT,flow(k+1)+λThighζT,fhigh(k+1)+λWlowζW,flow(k+1)+λWhighζW,fhigh(k+1)]$
where the λs are penalty parameters. The total supply airflow rate msa(k) used in Pfan(k) and Pcc(k) is given by $msa(k)=∑f∈Fmsa,f(k)=∑f∈F∑i∈Ifmsa,i(k)$. The optimal control commands are obtained by solving the minimization problem (10a) subject to constraints (10b) to (10r), which are described next.
The constraints due to the aggregate thermal dynamics of floor/meta-zone f are as follows:
$Tz,f(k+1)=Tz,f(k)+Δt[Toa(k)−Tz,f(k)τza,f+Tw,f(k)−Tz,f(k)τzw,f+Az,fηsol(k)+qint,f(k)+qac,f(k)Cz,f]$
(10b)
$Tw,f(k+1)=Tw,f(k)+Δt[Toa(k)−Tw,f(k)τwa,f+Tz,f(k)−Tw,f(k)τwz,f+Aw,fηsol(k)]$
(10c)
This is a discretized form of an RC (resistor-capacitor) network model, specifically a 2R2C model. The two states of the model are aggregate zone temperature (Tz,f) and wall temperature (Tw,f, a fictitious state). In constraint (10b), qac,f(k) is the heat influx due to the HVAC system and is given by $qac,f(k):=msa,f(k)Cpa(Tsa,f(k)−Tz,f(k))$. The model has seven parameters {Cz,f, τzw,f, τza,f, Az,f, τwz,f, τwa,f, Aw,f}. In the evaluation study, they are estimated using the algorithm presented in Ref. [20] and will be discussed later in Sec. 5.
The constraint due to the aggregate humidity dynamics of floor/meta-zone f is as follows:
$Wz,f(k+1)=Wz,f(k)+ΔtRgTz,f(k)VfPda×[ωint,f(k)+msa,f(k)Wca(k)−Wz,f(k)1+Wca(k)]$
(10d)
where Wz,f is the aggregate zone humidity ratio, Vf is the volume of meta-zone f, Rg is the specific gas constant of dry air, Pda is the partial pressure of dry air, and Wca is the conditioned air humidity ratio [33].
The constraints due to the control-oriented cooling and dehumidifying coil model, which was developed in our prior work [5], are as follows:
$Tca(k)=Tma(k)+mw(k)f(Tma(k),Wma(k),msa(k),mw(k))$
(10e)
$Wca(k)=Wma(k)+mw(k)g(Tma(k),Wma(k),msa(k),mw(k))$
(10f)
This specific functional form is chosen so that when the chilled water flowrate is zero, no cooling or dehumidification of the air occurs, so the conditioned air temperature and humidity ratio are equal to the mixed air temperature and humidity ratio: Tca = Tma and Wca = Wma, when mw = 0. The interested readers are referred to Ref. [5, Section 3.1.1] for details regarding the model.
The box constraints to maintain the temperature and humidity of the meta-zones within the allowed comfort limits are as follows:
$Tz,flow(k)−ζT,flow(k)≤Tz,f(k)≤Tz,fhigh(k)+ζT,fhigh(k)$
(10g)
$alowTz,f(k)+blow−ζW,flow(k)≤Wz,f(k)≤ahighTz,f(k)+bhigh+ζW,fhigh(k)$
(10h)
$ζT,flow(k+1),ζT,fhigh(k+1),ζW,flow(k+1),ζW,fhigh(k+1)≥0$
(10i)
The constraints (10g) and (10h) are softened using slack variables $ζT,flow(k)$, $ζT,fhigh(k)$, $ζW,flow(k)$, and $ζW,fhigh(k)$; constraint (10i) ensures that these slack variables are nonnegative. Imposing constraints directly on the relative humidity of zones (RHz) is difficult, as relative humidity is a highly nonlinear function of dry bulb temperature and humidity ratio [32, Chapter 1]. So, we linearize this function which gives us alow, blow, ahigh, and bhigh in Eq. (10h) and help in converting the constraints on relative humidity to humidity ratio (Wz).
The following constraint is for the outdoor airflow rate:
$moamin≤moa(k)≤moamax$
(10j)
where the minimum allowed value ($moamin$) is computed based on the ventilation requirements specified in ASHRAE 62.1 [34] and to maintain positive building pressurization.
In addition, there are several actuator constraints:
$Tca(k+1)≤min(Tca(k)+TcarateΔt,Tma(k+1),Tcahigh)$
(10k)
$Tca(k+1)≥max(Tca(k)−TcarateΔt,Tcalow)$
(10l)
$Wca(k)≤Wma(k)$
(10m)
$roa(k)=moa(k)/msa(k)$
(10n)
$roa(k+1)≤min(roa(k)+roarateΔt,roahigh)$
(10o)
$roa(k+1)≥max(roa(k)−roarateΔt,roalow)$
(10p)
$msa,flow≤msa,f(k)≤msa,fhigh$
(10q)
$Tca(k)≤Tsa,f(k)≤Tsahigh$
(10r)
Constraints (10k)(10l) and (10o)(10p) are to take into account the capabilities of the cooling coil and damper actuators, respectively. In constraints (10k) and (10m), the inequalities Tca(k + 1) ≤ Tma(k + 1) and Wca(k) ≤ Wma(k) ensure that the cooling coil can only cool and dehumidify the mixed air stream; it cannot add heat or moisture. Similarly, in constraint (10r) the inequality Tsa,f(k) ≥ Tca(k) ensures that the reheat coils can only add heat; they cannot cool. Constraint (10q) is to take into account the capabilities of the fan and aggregate capabilities of the VAV boxes. The limits $msa,flow$ and $msa,fhigh$ are computed using the VAV schedule from the mechanical drawings of a building as follows: $msa,flow:=∑i∈Ifmsa,ilow$ and $msa,fhigh:=∑i∈Ifmsa,ihigh$.

Overall, constraints (10b)(10d), (10g)(10i), and (10q)(10r) are $∀f∈F$. Constraints (10b)(10f), (10i)(10j), (10m)(10n), and (10q)(10r) are for k ∈ {j, …, j + NM − 1}, constraints (10g) and (10h) are for k ∈ {j + 1, …, j + NM}, and constraints (10k)(10l) and (10o)(10p) are for k ∈ {j − 1, …, j + NM − 2}. The control commands remain the same for M timesteps, as the control interval ΔT = MΔt, i.e., $uHLC(k)=uHLC(k+1)$ = … = $uHLC(k+M−1)$, $∀k∈{j,j+M,…,j+NM−1}$.

Note that of the states x(k) defined in Eq. (7), Tw,f is a fictitious state that cannot be measured, while the other two states aggregate zone temperature (Tz,f) and aggregate zone humidity ratio (Wz,f) are measured. So, we estimate the current state $x^(k)$ using a Kalman filter.

Comment 1. How best to aggregate zones into meta-zones is an open question. Too many meta-zones negate the benefit from aggregation. Too few will likely cause poor control performance since the difference in thermal properties and internal disturbances among various regions of the building will be lost when such a severe aggregation is performed. Works on model reduction of multi-zone buildings such as Ref. [35] can offer guidelines here. In our case study, we aggregate each floor of a building into a meta-zone as we expect some mass transfer among zones within a floor because of open doors, but not as much among floors.

### 3.2 Projection-Based Low-Level Controller.

The role of the LLC is to appropriately distribute the aggregate quantities—such as the total supply airflow rate and reheat power consumption—computed by the HLC to individual zones/VAVs. The LLC needs to do so by capturing two important properties: (i) it should consider the needs of individual zones and distribute accordingly and (ii) it should act in coherence with the HLC, so that there is minimal mismatch for the MPC optimization in the next round.

The LLC is a projection-based feedback controller that decides on the supply airflow rate and supply air temperature for each VAV box/zone. That is, the control command vector that the LLC needs to decide is:
$uLLC(k):=[msa,i(k),Tsa,i(k)]T∈ℜnz+nzrh$
where for msa,i, $i∈If,∀f∈F$, and for Tsa,i, $i∈Irh,f,∀f∈F$. It decides these control commands by using the following information from the HLC: (i) total allowed supply airflow rate to all the zones $msa(k)=∑f∈Fmsa,f(k)$, (ii) total allowed reheat power consumption $Preheat(k)=∑f∈FPreheat,f(k)$, (iii) the temperature at which the zones in each meta-zone should be maintained at Tz,f(k + 1), and (iv) the conditioned air temperature Tca(k). Here on in this section, we will be using the superscript HLC ($∙HLC$) for these variables to make it clear that these are obtained from the high-level controller.

First, the needs of each zone are assessed based on the current measured temperature Tz,i(k) and the range it should be in $[Tz,ihtg(k),Tz,iclg(k)]$ and are translated into the desired supply airflow rate $msa,id(k)$ and supply air temperature $Tsa,id(k)$. Then, these desired values along with the information obtained from the HLC are used to solve a projection problem to compute the control commands for all the zones, uLLC(k).

The procedure used to compute the desired values $msa,id(k)$ and $Tsa,id(k)$ is explained below. This is similar to the Dual Maximum control logic presented in Sec. 4; a schematic representation of it is shown in Fig. 6.

1. First the temperature range $[Tz,ihtg(k),Tz,iclg(k)]$ in which each zone should be is computed as follows: $Tz,ihtg(k)=max(Tz,fHLC(k+1)−T~zdb/2,Tz,flow)∀i∈If$ and $Tz,iclg(k)=min$$(Tz,fHLC(k+1)+T~zdb/2,Tz,fhigh)∀i∈If$, where $Tz,fHLC(k+1)$ is obtained from the HLC, $T~zdb$ is a deadband, and $Tz,flow$ and $Tz,fhigh$ are the limits used in constraint (10g).

2. If the zone temperature is between the cooling and heating setpoints ($Tz,i(k)∈[Tz,ihtg(k),Tz,iclg(k)]$), then the controller is in deadband mode. The supply airflow rate is desired to be at its minimum and no heating is required, i.e., $msa,id(k)=msa,ilow$ and $Tsa,id(k)=TcaHLC(k)$.

3. If the zone temperature is warmer than the cooling setpoint ($Tz,i(k)>Tz,iclg(k$)), then the controller is in cooling mode. The supply airflow rate is desired to be increased as needed and no heating is required, i.e., $msa,id(k)=min(msa,ilow+Km,iclg(Tz,i(k)−Tz,iclg(k)),msa,ihigh)$ and $Tsa,id(k)=TcaHLC(k)$.

4. If the zone temperature is cooler than the heating setpoint ($Tz,i(k))), then the controller is in heating mode. Heating is required, and the supply airflow rate is desired to be increased only if additional heating is needed, i.e., $Tsa,id(k)=min(TcaHLC(k)+KT,ihtg(Tz,ihtg(k)−Tz,i(k)),Tsahigh)$; if $Tsa,id(k)=Tsahigh$, then $msa,id(k)=min(msa,ilow+Km,ihtg(Tz,ihtg(k)$$−Tz,i(k)),msa,ihigh,reheat)$; otherwise, $msa,id(k)=msa,ilow$.

5. Finally, we impose the following rate constraints:
$msa,i(k−M)−msa,irateΔT≤msa,id(k)≤msa,i(k−M)+msa,irateΔTTsa,i(k−M)−Tsa,irateΔT≤Tsa,id(k)≤Tsa,i(k−M)+Tsa,irateΔT$
where msa,i(kM) and Tsa,i(kM) are the supply airflow rate and supply air temperature from the previous control time-step.

Fig. 6
Fig. 6
Close modal
These desired values—$msa,id(k)$ and $Tsa,id(k)$—along with information from the HLC are used to solve the following projection problem to obtain the control commands for all the zones, uLLC(k):
$minuLLC(k)∑f∈F∑i∈Ifλm,i(msa,i(k)−msa,id(k))2+∑f∈F∑i∈Irh,fλT,i(Tsa,i(k)−Tsa,id(k))2$
(11a)
subject to the following constraints:
$∑f∈F∑i∈Ifmsa,i(k)≤msaHLC(k)$
(11b)
$∑f∈F∑i∈Irh,fmsa,i(k)Cpa(Tsa,i(k)−TcaHLC(k))ηreheatCOPh≤PreheatHLC(k)$
(11c)
$msa,ilow≤msa,i(k)≤msa,ihigh,∀i∈If,∀f∈F$
(11d)
$TcaHLC(k)≤Tsa,i(k)≤Tsahigh,∀i∈Irh,f,∀f∈F$
(11e)
where the sets If and Irh,f are defined at the beginning of this section, λs are weights, $msaHLC(k)=∑f∈Fmsa,fHLC(k)$, and $PreheatHLC(k)=$$∑f∈FPreheat,fHLC(k)$.

Constraints (11b) and (11c) are to ensure that the total supply airflow rate and reheat power consumption do not exceed the limits computed by the HLC. Constraints (11d) and (11e) are to take in to account the capabilities of the VAV boxes and reheat coils. In constraint (11e), the inequality Tsa,i(k) ≥ Tca(k) ensures that the reheat coils can only add heat to the conditioned air and cannot cool. The upper limit on supply air temperature ($Tsahigh$) in constraint (11e) is to prevent thermal stratification [22].

## 4 Baseline Control (BL)

For zone climate control, we consider the Dual Maximum algorithm [22] as the baseline; a schematic representation of this algorithm is shown in Fig. 6. Even though Single Maximum is more commonly used, including in the Innovation Hub building, we choose Dual Maximum for the baseline, as it is more energy-efficient among the two [22,36]. The Dual Maximum controller operates in three modes based on the zone temperature (Tz,i): (i) cooling, (ii) heating, and (iii) deadband. The zone’s supply airflow rate (msa,i) and supply air temperature (Tsa,i) are varied based on the mode as explained below.

1. Cooling mode: If the zone temperature is warmer than the cooling setpoint, then the controller is in cooling mode. The supply airflow rate (msa,i) is varied between the minimum ($msa,ilow$) and maximum ($msa,ihigh$) as needed, and the supply air temperature (Tsa,i) is equal to the conditioned air temperature (Tca), i.e., no reheat.

2. Heating mode: If the zone temperature is below the heating setpoint, then the controller is in heating mode. First, the supply air temperature (Tsa,i) is increased up to the maximum allowed value ($Tsahigh$) as needed to maintain the zone temperature at the heating setpoint. If the zone temperature still cannot be maintained at the heating setpoint, then the supply airflow rate is increased between the minimum ($msa,ilow$) and the heating maximum ($msa,ihigh,reheat$) values.

3. Deadband mode: If the zone temperature is between the heating and cooling setpoints, then the controller is in deadband mode. The supply airflow rate is kept at the minimum, and the supply air temperature is equal to the conditioned air temperature, i.e., no reheat.

In the case of VAV boxes that do not have reheat coils, the logic during cooling and deadband modes are the same. In heating mode, the supply airflow rate is at the minimum and the supply air temperature is equal to the conditioned air temperature, as the VAV box cannot heat.

For the AHU-level commands, the BL controller uses fixed conditioned air temperature that is determined based on expected thermal (sensible and latent) load, and fixed outdoor airflow rates based on ventilation requirements, e.g., ASHRAE 62.1 [34]. Another consideration in choosing outdoor airflow rate is building positive pressurization requirements [32].

## 5 Simulation Setup

Recall that the plant is based on an air handling unit serving 33 zones, of which 29 are equipped with reheat coils, and the remaining four do not have reheat coils (cooling only). Of the 29 VAV boxes with reheat, three of them serve laboratories which are equipped with fume hoods (209, 303, and 310), and one of them serve restrooms (103). The VAV boxes serving these labs need to be controlled to satisfy the negative pressurization requirements with respect to corridor, so we assume that they operate according to the existing rule based feedback control strategy. Therefore, nz = 29 and $nzrh=25$; for msa,i, $i∈{I1:={101–102,104–109$}, $I2:={201–208,210–212}$, $I3:={301–302,304–309$, 311–312} and for Tsa,i, $i∈{Irh,1:={I1∖{106,107}}$, $Irh,2:={I2$$∖{205}}$, $Irh,3:={I3∖{307}}}$. The sets I1, I2, and I3 defined above are the VAVs/zones in floors 1, 2, and 3, respectively. The sets Irh,1, Irh,2, and Irh,3 exclude the VAV boxes which do not have a reheat coil.

The outdoor weather data used in simulations is obtained from the National Solar Radiation Database2 for Gainesville, FL. As mentioned in Sec. 2.1, the internal heat load due to occupants are computed based on the number of occupants provided to the zone block. We assume that the zones are occupied from Monday to Friday between 8:00 a.m. to noon and 1:00 p.m. to 5:00 p.m., with the total number of occupants (np,f) in floor 1 as 24, in floor 2 as 26, and in floor 3 as 22. We assume a power density of 12.92 W/m2 (1.2 W/ft2) for internal heat load due to lighting and equipment. For special purpose rooms like electrical and telecommunication, we use a higher power density of 53.82 W/m2 (5 W/ft2). These heat loads from lighting and equipment are assumed to be halved during weekends.

The following zone temperature and humidity limits are used in the simulations: $Tzlow=21.1∘C(70∘F)$, $Tzhigh=23.3∘C(74∘F)$, $RHzlow=20%$, and $RHzhigh=65%$. The chosen thermal comfort envelope is shown in Fig. 7. Typically, the zone temperature limits during unoccupied mode (unocc) are relaxed when compared to the occupied mode (occ), i.e., $[Tzlow,occ,$$Tzhigh,occ]⊆[Tzlow,unocc,Tzhigh,unocc]$. Due to its usage, the Innovation Hub building is always operated in occupied mode, so we assume the same in simulations. For the simulation results reported later, the zone temperature violation is computed as $max(Tz(k)−Tzhigh,Tzlow−Tz(k),0)$ and the zone relative humidity violation is computed as $max(RHz(k)−RHzhigh,RHzlow−RHz(k),0)$, with the upper and lower limits mentioned above.

Fig. 7
Fig. 7
Close modal

The fan power coefficient αfan in Eq. (2) is 14.2005 W/(kg/s)3, which is obtained using a least squares fit to data collected from the building. The parameters of the cooling and dehumidifying coil model used in the plant are fit using the procedure explained in Sec. 2.1.2 of Ref. [5]. The root-mean-square errors for the validation data set are 0.25 °C (0.46 °F, 2%) for Tca and 0.22 × 10−4 kgw/kgda (2.6%) for Wca.

MZHC parameters: The optimization problems in HLC and LLC are solved using CasADi [37] and IPOPT [38], a nonlinear programming (NLP) solver, on a Desktop Windows computer with 16 GB RAM and a 3.60 GHZ × 8 CPU. As mentioned in Sec. 3.1, Δt = 5 min, ΔT = 15 min, T = 24 h, N = 96, and M = 3. The number of decision variables for the HLC are 7008 and for the LLC are 54. On an average, it takes only 3.28 s to solve the optimization problem in the HLC and 0.018 s to solve the optimization problem in the LLC.

The parameters for the control-oriented cooling and dehumidifying coil model are fit using the procedure explained in Ref. [5]. For the validation data set, the root-mean-square errors are 0.97 °C (1.75 °F, 7.6%) for Tca and 0.63 × 10−4 kgw/kgda (7.6%) for Wca.

Since the Innovation Hub building has three floors, we aggregate it into three meta-zones, i.e., $f∈{1,2,3}=:F$. The parameters of the aggregate thermal dynamics model for each meta-zone are estimated using the algorithm presented in Ref. [20]. The interested reader is referred to Ref. [30] for the values of the parameters and out-of-sample prediction accuracy of the identified model.

For the aggregate humidity dynamics model, floor volumes used are V1 = 1036.6 m3, V2 = 1504.1 m3, and V3 = 1330.8 m3, which are obtained from mechanical drawings of the building.

The following limits are used for the zone temperature constraint (10g): $Tz,flow=$21.1 °C (70 °F) and $Tz,fhigh=$23.3 °C (74 °F). The coefficients for the humidity constraint in (10h) are ahigh = 0.000621 kgw/kgda/ °C and bhigh = −0.173323 kgw/kgda, which corresponds to a relative humidity of 60%, and alow = 0.000203 kgw/kgda/ °C and blow = −0.056516 kgw/kgda, which corresponds to a relative humidity of 20%. We introduce a factor of safety by using a slightly tighter higher limit of 60% for the relative humidity of the zones when compared to the thermal comfort envelope presented in Fig. 7.

The minimum allowed value for the outdoor airflow rate ($moamin$) is 3.24 kg/s (5700 cfm), which is obtained from the AHU schedule in the mechanical drawings for the building. The maximum possible value for the outdoor airflow rate ($moamax$) is 8.52 kg/s (15 000 cfm). The various limits on the supply airflow rates are obtained using the VAV schedule for the building. The remaining limits used in the controllers are as follows: $roalow=0%$, $roahigh=100%$, $T~zdb=0.56∘C$ (1 °F), $Tcalow=11.67∘C$ (53 °F), $Tcahigh=17.2∘C$ (63 °F), and $Tsahigh=30∘C$ (86 °F). The higher limit on the conditioned air temperature ($Tcahigh$) is to introduce a factor of safety and make the controller robust.

The MPC controller requires predictions of the various exogenous inputs specified in Eq. (9). We compute the loads due to occupants in qint,f and ωint,f based on the occupancy profile used in simulating the plant. The outdoor weather related exogenous inputs are assumed to be fully known.

BL parameters: The cooling and heating setpoints are chosen to be 21.1 °C (70 °F) and 23.3 °C (74 °F), respectively. The minimum, maximum, and heating maximum values for the supply airflow rate of all the VAV boxes are obtained from the VAV schedule for the building. The maximum allowed value for the supply air temperature ($Tsahigh$) is 30 °C (86 °F). The conditioned air temperature (Tca) is kept at a constant value of 11.67 °C (53 °F). Typically, Tca is kept at 12.8 °C (55 °F), especially in hot-humid climates, to ensure humidity control [39], but we assume there is a 1.11 °C (2 °F) increase in temperature because of the heat from the draw-through supply fan in the AHU, so we keep it at 11.67 °C (53 °F) to compensate for it. The outdoor airflow rate is kept at 3.24 kg/s (5700 cfm), which is obtained from mechanical drawings for the building.

## 6 Results and Discussion

Performance of the controllers is compared using three types of outdoor weather conditions: hot-humid (Jul/06 to Jul/13), mild (Feb/19 to Feb/26), and cold (Jan/30 to Feb/06). The proposed controller reduces energy use significantly compared to BL, from approximately 11% to 68% depending on weather, see Fig. 8. The indoor climate control performance of MZHC and BL is similar. With MZHC, the RMSE (root-mean-square error) of zone temperature violation is 0.1 °C (0.18 °F) and the RMSE of zone relative humidity violation is 0.05%, while with BL they are 0.01 °C (0.02 °F) and 0%, respectively. The computational cost of the proposed MZHC is small, just a few seconds are needed to compute decisions at every control update. On an average it takes 3.28 s to solve the optimization problem in the HLC and 0.018 s to solve the optimization problem in the LLC.

Fig. 8
Fig. 8
Close modal

Simulation results for the different weather conditions are discussed in detail next.

### 6.1 Hot-Humid Week.

Figure 9 shows the simulation time traces for a hot-humid week. Both controllers lead to negligible violations of aggregate zone temperature (Tz,f) and relative humidity (RHz,f) constraints. Data for the three meta-zones are shown in Fig. 9(c), and for three individual zones, one from each floor, are shown in Fig. 10. BL ensures that dry enough air is supplied to the zones at all times by keeping the conditioned air temperature (Tca) at a constant low value of 11.67 °C (53 °F), and hence, the humidity limit is not violated. In the case of MZHC, the humidity constraint is found to be active always. This can be seen in Fig. 9(c); recall that we use a tighter constraint of 60% instead of 65% to introduce a factor of safety. This active constraint limits the increase in Tca, which can be seen in Fig. 9(d). One of the reasons reheating is required even during such a hot week (the outdoor air temperature is as high as 32.2 °C/90 °F) is because of this active humidity constraint, which requires dry, and thus, cold air to be supplied to the zones. This could also be one of the reasons MZHC decides to maintain the zone temperatures at the lower limit (see Fig. 9(c)). Keeping the zones at a higher temperature will require an increase in reheating energy.

Fig. 9
Fig. 9
Close modal
Fig. 10
Fig. 10
Close modal

Most of the prior works on using MPC for HVAC control ignore humidity and latent heat in their formulation. In an attempt to reduce energy/cost, such controllers are likely to make decisions during these hot-humid conditions which will lead to poor indoor humidity, as they are unaware of the factors mentioned above [5]. Thus, such controllers cannot be used particularly in hot-humid climates.

The energy savings by MZHC is because of two main reasons. One, MZHC increases the Tca as long as the humidity constraints are not violated, while BL uses a conservatively designed value as explained above. This leads to a reduction in cooling energy consumption by MZHC, see Figs. 8 and 9(b)). Two, the warmer Tca supplied to the VAV boxes requires lesser reheating in the case of MZHC. This leads to a reduction in the reheat energy consumption, which can be seen in Figs. 8 and 9(b). The decisions regarding the supply airflow rates are found to be the same for both the controllers, and thus, the fan energy consumptions are identical.

Since the outdoor air is hot and humid (see Fig. 9(a)), bringing in more than the minimum outdoor airflow rate (moa) required will increase the sensible and latent loads on the cooling coil. So, MZHC decides to keep moa at the minimum, see Fig. 9(d).

### 6.2 Mild and Cold Weeks.

As seen from Fig. 8, energy savings are higher than 60% during mild and cold weeks. This significant reduction in energy consumption can be attributed to three main reasons. Two of them are the same as those explained in Sec. 6.1, but, the effects here are more prominent. For instance, unlike the hot-humid week, the humidity constraint is only intermittently active since the outdoor weather is drier. This provides more room for optimizing the conditioned air temperature, which has two important implications: (i) reduction in the cooling energy consumption and (ii) minimal reheat energy consumption. The third is that, since the outdoor weather is cold and dry, MZHC also decides to use “free” cooling when possible by bringing in more than the minimum outdoor air required which leads to further reduction in cooling energy consumption.

We omit the detailed results for the mild and cold weeks here; the interested reader can find them in Ref. [30].

Comment 2. Recall that as mentioned in Sec. 4, the Innovation Hub building uses Single Maximum algorithm for zone climate control as opposed to the Dual Maximum algorithm presented here [22]. In the case of Single Maximum algorithm, the minimum allowed value for the supply airflow rate has to be high enough so that the heating load can be met with low enough supply air temperature to prevent thermal stratification [22]. While in the case of Dual Maximum algorithm, the airflow rate is varied during the heating mode; thus, the minimum allowed airflow rate is not limited by stratification. Therefore, using Single Maximum algorithm leads to higher fan, cooling, and reheating energy consumptions. This also implies that the energy savings by the proposed controller will be even higher.

## 7 Conclusion

The proposed control architecture is designed to address a number of limitations in the existing literature on multi-zone building control using MPC. The main one is the reliance on a high-resolution multi-zone model, which can be challenging to obtain. A low-resolution model of the building is more convenient since such a model can be identified in a tractable manner from measurements. The challenge then is to convert the MPC decisions, which are computed for the fictitious zones in the model, to the commands for the VAV boxes of the actual building. The proposed architecture does that by posing this conversion as a projection problem that uses not only what the MPC computes but also feedback from zones. The result is a principled method of computing VAV box commands that are consistent with the optimal decisions made by the MPC without needing dynamic models of individual zones. At the same time, the use of feedback from the zones ensures that zone climate states are also close to the aggregate climate states computed by the MPC. Additionally, our formulation takes into account humidity and latent heat, something that is often ignored in the literature on MPC for HVAC systems.

The positive simulation results provide confidence on the effectiveness of the proposed controller especially because of the large plant model mismatch. The simulation testbed mimics a real building closely, including the heterogeneous nature of the zones in the building. Another observation from the simulations that should be emphasized is the need to incorporate humidity and latent heat. The indoor humidity constraint is seen to be active when using the proposed controller, especially during hot-humid weather. Without humidity being explicitly considered, the controller is likely to have caused high space humidity in an effort to reduce energy use.

There are many avenues for further exploration, such as experimental verification, modifying the formulation to minimize cost instead of energy and to include demand charges, improving methods to forecast disturbances, etc.

## Acknowledgment

This research reported here has been partially supported by the NSF through award # 1934322 and the State of Florida through a REET (Renewable Energy and Energy Efficient Technologies) grant. We thank David Brooks for help in understanding the design and operation of Innovation Hub’s HVAC system.

## Conflict of Interest

There are no conflicts of interest.

## Data Availability Statement

Data provided by a third party listed in Acknowledgment.

## References

1.
Serale
,
G.
,
Fiorentini
,
M.
,
Capozzoli
,
A.
,
Bernardini
,
D.
, and
,
A.
,
2018
, “
Model Predictive Control (MPC) for Enhancing Building and HVAC System Energy Efficiency: Problem Formulation, Applications and Opportunities
,”
Energies
,
11
(
3
), p.
631
.
2.
Shaikh
,
P. H.
,
Nor
,
N. B. M.
,
Nallagownden
,
P.
,
Elamvazuthi
,
I.
, and
Ibrahim
,
T.
,
2014
, “
A Review on Optimized Control Systems for Building Energy and Comfort Management of Smart Sustainable Buildings
,”
Renewable. Sustainable. Energy. Rev.
,
34
, pp.
409
429
.
3.
Goyal
,
S.
, and
Barooah
,
P.
,
2013
, “
Energy-Efficient Control of an Air Handling Unit for a Single-Zone VAV System
,”
IEEE Conference on Decision and Control
,
Florence, Italy
, pp.
4796
4801
.
4.
Joe
,
J.
, and
Karava
,
P.
,
2019
, “
A Model Predictive Control Strategy to Optimize the Performance of Radiant Floor Heating and Cooling Systems in Office Buildings
,”
Appl. Energy.
,
245
, pp.
65
77
.
5.
Raman
,
N. S.
,
,
K.
,
Chen
,
B.
,
Ingley
,
H. A.
, and
Barooah
,
P.
,
2020
, “
Model Predictive Control for Energy-Efficient HVAC Operation With Humidity and Latent Heat Considerations
,”
Appl. Energy.
,
279
(Dec.), p.
115765
.
6.
Chen
,
X.
,
Wang
,
Q.
, and
Srebric
,
J.
,
2016
, “
Occupant Feedback Based Model Predictive Control for Thermal Comfort and Energy Optimization: A Chamber Experimental Evaluation
,”
Appl. Energy.
,
164
, pp.
341
351
.
7.
Ma
,
J.
,
Qin
,
J.
,
Salsbury
,
T.
, and
Xu
,
P.
,
2012
, “
Demand Reduction in Building Energy Systems Based on Economic Model Predictive Control
,”
Chem. Eng. Sci.
,
67
(
1
), pp.
92
100
. Dynamics, Control and Optimization of Energy Systems.
8.
Bengea
,
S. C.
,
Kelman
,
A. D.
,
Borrelli
,
F.
,
Taylor
,
R.
, and
Narayanan
,
S.
,
2013
, “
Implementation of Model Predictive Control for An HVAC System in a Mid-size Commercial Building
,”
HVAC &R Res.
,
20
(
1
), pp.
121
135
.
9.
,
N.
,
Su
,
Y.
,
Su
,
R.
, and
Poolla
,
K.
,
2016
, “
Token Based Scheduling for Energy Management in Building HVAC Systems
,”
Appl. Energy.
,
173
, pp.
67
79
.
10.
Png
,
E.
,
Srinivasan
,
S.
,
Bekiroglu
,
K.
,
Chaoyang
,
J.
,
Su
,
R.
, and
Poolla
,
K.
,
2019
, “
An Internet of Things Upgrade for Smart and Scalable Heating, Ventilation and Air-conditioning Control in Commercial Buildings
,”
Appl. Energy.
,
239
, pp.
408
424
.
11.
Ma
,
Y.
,
Richter
,
S.
, and
Borrelli
,
F.
,
2012
, “
Chapter 14: Distributed Model Predictive Control for Building Temperature Regulation
,”
Control Optim. Differ. Algebraic Constraints
,
22
(March), pp.
293
314
.
12.
Mei
,
J.
, and
Xia
,
X.
,
2018
, “
Multi-Zone Building Temperature Control and Energy Efficiency Using Autonomous Hierarchical Control Strategy
,”
2018 IEEE 14th International Conference on Control and Automation (ICCA)
,
Anchorage, AL
, pp.
884
889
.
13.
Patel
,
N. R.
,
Risbeck
,
M. J.
,
Rawlings
,
J. B.
,
Wenzel
,
M. J.
, and
Turney
,
R. D.
,
2016
, “
Distributed Economic Model Predictive Control for Large-scale Building Temperature Regulation
,”
2016 American Control Conference (ACC)
,
Boston, MA
, pp.
895
900
.
14.
Yang
,
Y.
,
Hu
,
G.
, and
Spanos
,
C. J.
,
2020
, “
HVAC Energy Cost Optimization for a Multizone Building Via a Decentralized Approach
,”
IEEE Transactions on Automation Science and Engineering
,
17
(
4
), pp.
1950
1960
.
15.
Long
,
Y.
,
Liu
,
S.
,
Xie
,
L.
, and
Johansson
,
K. H.
,
2016
, “
A Hierarchical Distributed MPC for HVAC Systems
,”
2016 American Control Conference (ACC)
,
Boston, MA
, pp.
2385
2390
.
16.
Jorissen
,
F.
, and
Helsen
,
L.
,
2016
, “
Towards An Automated Tool Chain for MPC in Multi-Zone Buildings
,”
4th International Conference on High Performance Buildings
,
West Lafayette, IN
, pp.
1
10
.
17.
Li
,
X.
, and
Wen
,
J.
,
2014
, “
Review of Building Energy Modeling for Control and Operation
,”
Renew. Sustain. Energy Rev.
,
37
, pp.
517
537
.
18.
Kim
,
D.
,
Cai
,
J.
,
Ariyur
,
K. B.
, and
Braun
,
J. E.
,
2016
, “
System Identification for Building Thermal Systems Under the Presence of Unmeasured Disturbances in Closed Loop Operation: Lumped Disturbance Modeling Approach
,”
Build. Environ.
,
107
, pp.
169
180
.
19.
Zeng
,
T.
, and
Barooah
,
P.
,
2020
, “
Identification of Network Dynamics and Disturbance for a Multi-zone Building
,”
IEEE Trans. Control Syst. Technol.
,
28
(Aug.), pp.
2061
2068
.
20.
Guo
,
Z.
,
Coffman
,
A. R.
,
Munk
,
J.
,
Im
,
P.
,
Kuruganti
,
T.
, and
Barooah
,
P.
,
2020
, “
Aggregation and Data Driven Identification of Building Thermal Dynamic Model and Unmeasured Disturbance
,”
Energy Build
,
231
, p.
110500
.
21.
Fritzson
,
P.
, and
Engelson
,
V.
,
1998
,
Modelica — A Unified Object-Oriented Language for System Modeling and Simulation
,
E.
Jul
, ed.,
Springer
,
Berlin/Heidelberg, Germany
, pp.
67
90
.
22.
ASHRAE
,
2011
,
The ASHRAE Handbook: HVAC Applications
, SI ed.,
ASHRAE
,
Atlanta, GA
.
23.
Jana
,
R. L.
,
Dey
,
S.
, and
Dasgupta
,
P.
,
2020
, “
A Hierarchical HVAC Control Scheme for Energy-Aware Smart Building Automation
,”
Assoc. Comput. Mach.
,
25
(
4
), pp.
31:1
31:33
.
24.
Liang
,
W.
,
Quinte
,
R.
,
Jia
,
X.
, and
Sun
,
J.-Q.
,
2015
, “
MPC Control for Improving Energy Efficiency of a Building Air Handler for Multi-Zone VAVs
,”
Build. Environ.
,
92
, pp.
256
268
.
25.
Ma
,
Y.
,
Matuško
,
J.
, and
Borrelli
,
F.
,
2015
, “
Stochastic Model Predictive Control for Building HVAC Systems: Complexity and Conservatism
,”
IEEE Trans. Control Syst. Technol.
,
23
(
1
), pp.
101
116
.
26.
Huang
,
S.
,
Lin
,
Y.
,
Chinde
,
V.
,
Ma
,
X.
, and
Lian
,
J.
,
2021
, “
Simulation-Based Performance Evaluation of Model Predictive Control for Building Energy Systems
,”
Appl. Energy.
,
281
, p.
116027
.
27.
Jorissen
,
F.
,
Reynders
,
G.
,
Baetens
,
R.
,
Picard
,
D.
,
Saelens
,
D.
, and
Helsen
,
L.
,
2018
, “
Implementation and Verification of the IDEAS Building Energy Simulation Library
,”
J. Build. Perform. Simulat.
,
11
(
6
), pp.
669
688
.
28.
IBPSA
, “Modelica IBPSA library,” https://github.com/ibpsa/modelica-ibpsa.
29.
Wetter
,
M.
,
Zuo
,
W.
,
Nouidui
,
T. S.
, and
Pang
,
X.
,
2014
, “
Modelica Buildings Library
,”
J. Build. Perform. Simulat.
,
7
(
4
), pp.
253
270
.
30.
Raman
,
N. S.
,
Chaturvedi
,
R. U.
,
Guo
,
Z.
, and
Barooah
,
P.
,
2021
, “
MPC-Based Hierarchical Control of a Multi-Zone Commercial HVAC System
,” arXiv preprint:2102.02914.
31.
Roulet
,
C.-A.
,
Heidt
,
F.
,
,
F.
, and
Pibiri
,
M.-C.
,
2001
, “
Real Heat Recovery With Air Handling Units
,”
Energy Build.
,
33
(
5
), pp.
495
502
.
32.
ASHRAE
,
2017
,
The ASHRAE Handbook Fundamentals
, SI ed.,
ASHRAE
,
Atlanta, GA
.
33.
Goyal
,
S.
, and
Barooah
,
P.
,
2012
, “
A Method for Model-Reduction of Non-Linear Thermal Dynamics of Multi-Zone Buildings
,”
Energy Build.
,
47
(April), pp.
332
340
.
34.
ASHRAE
,
2016
, ANSI/ASHRAE standard 62.1-2016, ventilation for acceptable air quality.
35.
Deng
,
K.
,
Goyal
,
S.
,
Barooah
,
P.
, and
Mehta
,
P. G.
,
2014
, “
Structure-Preserving Model Reduction of Nonlinear Building Thermal Models
,”
Automatica
,
50
(
4
), pp.
1188
1195
.
36.
Goyal
,
S.
,
Ingley
,
H.
, and
Barooah
,
P.
,
2013
, “
Occupancy-based Zone Climate Control for Energy Efficient Buildings: Complexity Vs. Performance
,”
Appl. Energy.
,
106
, pp.
209
221
.
37.
,
J. A. E.
,
Gillis
,
J.
,
Horn
,
G.
,
Rawlings
,
J. B.
, and
Diehl
,
M.
,
2019
, “
Casadi: A Software Framework for Nonlinear Optimization and Optimal Control
,”
Mathe. Program. Comput.
,
11
(
1
), pp.
1
36
.
38.
Wächter
,
A.
, and
Biegler
,
L. T.
,
2006
, “
On the Implementation of An Interior-Point Filter Line-search Algorithm for Large-Scale Nonlinear Programming
,”
Math. Program.
,
106
(
1
), pp.
25
57
.
39.
Williams
,
J.
,
2013
, “
Why is the Supply Air Temperature 55F?
http://8760engineeringblog.blogspot.com/2013/02/why-is-supply-air-temperature-55f.html, Accessed August 3, 2020.