MPC-Based Hierarchical Control of a Multi-Zone Commercial HVAC System
Naren Srivaths Raman, Rahul Umashankar Chaturvedi, Zhong Guo, Prabir Barooah
MMPC-Based Hierarchical Control of a Multi-ZoneCommercial HVAC System
Naren Srivaths Raman ∗ , Rahul Umashankar Chaturvedi, Zhong Guo, and Prabir Barooah Department of Mechanical and Aerospace EngineeringUniversity of FloridaGainesville, FL 32611
This paper presents a novel architecture for model predictivecontrol (MPC) based indoor climate control of multi-zonebuildings to provide energy efficiency. Unlike prior workswe 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 buildingwhich is divided into a small number of “meta-zones” thatcan be easily identified using existing data-driven modelingtechniques. The proposed architecture is hierarchical. Atthe higher level, an MPC controller uses the low-resolutionmodel to make decisions for the air handling unit (AHU) andthe meta-zones. Since the meta-zones are fictitious, a lowerlevel controller converts the high-level MPC decisions intocommands for the individual zones by solving a projectionproblem that strikes a trade-off between two potentially con-flicting goals: the AHU-level decisions made by the MPC arerespected while the climate of the individual zones is main-tained within the comfort bounds. The performance of theproposed 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 mul-tiple weather conditions show the effectiveness of the pro-posed controller in terms of energy savings, climate control,and computational tractability.
The application of model predictive control (MPC)for commercial heating, ventilation, and air conditioning(HVAC) systems for both energy efficiency and demand flex-ibility has been an active area of research; see the review ar-ticles [1, 2] and the references therein.Several of the MPC formulations proposed in the pastare for buildings with one zone [3, 4, 5, 6] or a small num-ber of zones [7, 8]. A direct extension of such formula-tions for large multi-zone buildings has two main challenges.First, solving the underlying optimization problem in MPCfor a building with a large number of zones is computa-tionally complex because of the large number of decision ∗ Corresponding author. Email: narensraman@ufl.edu variables. To reduce the computational complexity, sev-eral distributed and hierarchical approaches have been pro-posed [9, 10, 11, 12, 13, 14, 15]. The second challenge, whichhas attracted far less attention, is that MPC requires a “high-resolution” model of the thermal dynamics of a multi-zonebuilding. High-resolution means that the temperature of ev-ery zone in the building is a state in the model and the controlcommands for every zone are inputs in the model. One wayof obtaining such a multi-zone model is by first construct-ing a “white box” model, such as by using a building energymodeling software, and then simplifying it to make it suit-able for MPC, e.g., [16]. But constructing a white box modelis expensive; it requires significant effort [17]. Moreover,the resulting model may not reflect the building as is. An-other way of obtaining a high-resolution multi-zone model isby utilizing data-driven techniques, which use input-outputmeasurements. Getting reliable estimates using data-drivenmodeling is challenging even for a single-zone building, as abuilding’s thermal dynamics is affected by a non-trivial andunmeasurable disturbance, the heat gain from occupants andtheir use of equipment, that strongly affects quality of theidentified model [18]. In the case of multi-zone model identi-fication, it becomes intractable since the model has too manydegrees of freedom: as many unknown disturbance signalsas there are number of zones. To the best of our knowledge,there are no works on reliable identification of multi-zonebuilding models without making assumptions on the natureof the disturbance affecting individual zones [19].In addition to the challenges mentioned above, most ofthe prior works—whether on single-zone or on multi-zonebuildings—ignore humidity and latent heat in their MPC for-mulations. The inclusion of moisture requires a computa-tionally convenient cooling and dehumidifying coil model.MPC formulations which exclude humidity can lead to poorhumidity control, or higher energy usage as they are unawareof the latent load on the cooling coil [5].In this work, we propose a humidity-aware MPC for-mulation for a multi-zone building with a variable air vol-ume (VAV) HVAC system. Figure 1 shows the schematic ofsuch a system.To overcome the challenges mentioned above, we pro-1 a r X i v : . [ ee ss . S Y ] F e b eturn airExhaustairOutdoor air Mixedair Cooling and dehumidifying coil
Tca
Chilledwater msaTma
Fan
Toa WcaWmaWoa Twimw Tsa,1mramoa Zone 2Zone 1
VAVW/REHEAT msa,1Tsa,2msa,2
VAVW/REHEAT
Tz,1 W z,1
Tz,2 W z,2 pduct
Air Handling Unit
TraWra Wsa,1Wsa,2
Fig. 1: Schematic of a multi-zone—specifically, a two zone—commercial variable-air-volume HVAC system. In this figure,oa: outdoor air, ra: return air, ma: mixed air, ca: conditioned air, and sa: supply air.pose a two-level control architecture. The high-level con-troller (HLC) decides on the AHU-level control commands.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 partof the building consisting of several zones. In the case studypresented here, a 33 zone three-floor building is aggregatedto a 3 meta-zone model, with each meta-zone correspondingto 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, inwhich disturbance in all the zones are aggregated into onesignal, can be directly identified from measurements col-lected from the building. The identification problem of sucha single-zone equivalent model is more tractable [20]. In thispaper, we use the system identification method from [20],though other identification methods can also be used. Sincethe HLC uses a low-resolution model with a much smallernumber of meta-zones than that in the building, its compu-tational complexity is low. However, this reduction of com-putational complexity creates a different challenge. Sincethe decision variables of the optimization problem in theHLC correspond to the meta-zones (air flow rate, temper-ature, etc.), they do not correspond to those for the actualzones of the building. The low-level controller (LLC) is nowused to compute the control commands for individual zones.It does so by solving a projection problem that appropriatelydistributes aggregate quantities computed by the HLC to in-dividual zones. The LLC uses feedback from each zone toassess their needs and ensures indoor climate of each zone ismaintained.The proposed controller—that includes the HLC andLLC–is hereafter referred to as
MZHC which stands for multi-zone hierarchical controller . Its performance is as-sessed through simulations on a “virtual building” plant. Theplant is representative of a section of the Innovation Hubbuilding comprising of 33 zones and is located at the Uni-versity of Florida campus. The plant is constructed usingModelica [21]. The performance of the proposed controlleris compared with that of
Dual Maximum controller as a base-line [22]. The dual maximum controller—which is referredto as BL (for baseline)—is a rule-based controller, and is one of the more energy efficient controllers among those used inpractice [22]. Simulation results show that using the pro-posed controller provides significant energy savings whencompared to BL while maintaining indoor climate.Compared to the literature on MPC design for multi-zone building HVAC systems, our work makes four principalcontributions, with details discussed in Section 1.1. (i) Thefirst contribution is that the proposed method does not as-sume availability of a high-resolution model of the multi-zone building which is difficult to obtain. Instead, it canutilize existing data-driven methods that can quickly iden-tify a low-resolution model of the multi-zone building frommeasurements. (ii) Since the MPC part of the proposed con-troller uses a low-resolution model with a small number ofmeta-zones, the method is scalable to buildings with a largenumber of zones. Although distributed iterative computationhas been proposed in the literature as an alternate approachto reducing computational complexity, ours can be solved ina centralized setting. (iii) The third contribution is the incor-poration of humidity and latent heat in our multi-zone MPCformulation, which has been largely ignored in the literatureon MPC for buildings, and especially so in the literature onmulti-zone building MPC. Our simulations show that whenusing MZHC , the indoor humidity constraint is active, es-pecially during hot-humid weather. Without humidity beingexplicitly considered, the controller would have caused highspace humidity in an effort to reduce energy use. (iv) Thefourth contribution is a realistic evaluation of the proposedcontroller in a high-fidelity simulation platform that intro-duces a large plant-model mismatch. In many prior workson multi-zone MPC, the model used by the controller is thesame as that used in simulating the plant. In contrast, theonly information provided to the proposed controller aboutthe building is sensor measurements (past data for modelidentification and real-time data during control) and designparameters such as expected occupancy, minimum designairflow rates for each VAV box, etc.The rest of this paper is organized as follows. Sec-tion 1.1 discusses our work in relation to the literature onmulti-zone MPC. Section 2 describes a multi-zone buildingequipped with a VAV HVAC system and the models we usein simulating the plant (the system to be controlled). Sec-2ion 3 presents the proposed MPC-based hierarchical con-troller. Section 4 describes a rule-based baseline controllerwith which the performance of the proposed controller iscompared. The simulation setup is described in Section 5.Simulation results are presented and discussed in Section 6.Finally, the main conclusions are provided in Section 7.
Several distributed and hierarchical approaches havebeen proposed to reduce the computational complexity ofMPC for multi-zone buildings [9, 10, 11, 12, 13, 14, 15].In [9], a hierarchical distributed algorithm called token-basedscheduling is proposed to vary the supply airflow rate to thezones. A modified version of this algorithm is used in [10] tominimize the energy consumption of a multi-zone buildinglocated at the Nanyang Technological University, Singaporecampus.In [12], a two-layered control architecture is proposedfor operating a VAV HVAC system. The upper layer is anopen loop controller, while the lower layer is based on MPCand it varies the supply airflow rates to the zones. Similarto [9], [10], and [12], the works [14, 13, 15] consider vary-ing only the zone-level control inputs such as the supply air-flow rates and zone temperature set points. These works ex-clude the AHU-level control inputs such as the conditionedair temperature and outside airflow rate. Unlike the worksmentioned above, the work [23] uses MPC to vary only theAHU-level control inputs; the zone-level control inputs areexcluded in this formulation.One of the few works similar to ours is [11], as theyconsider both the zone-level and AHU-level control inputsin their formulation. But their algorithm requires a high-resolution multi-zone model, and they do not consider hu-midity and latent heat in their formulation.
Our focus is a multi-zone building equipped with avariable-air-volume (VAV) HVAC system, whose schematicis shown in Figure 1. In such a system, part of the air ex-hausted from the zones is recirculated and mixed with out-door air. This mixed air is sent through the cooling coilwhere the air is cooled and dehumidified to the conditionedair temperature ( T ca ) and humidity ratio ( W ca ). This condi-tioned air is then sent through the supply ducts to the VAVboxes, which have a damper to control airflow, and finallysupplied to the zones. Some VAV boxes have reheat coils;they can change temperature of supply air but not humid-ity, i.e., T sa , i ≥ T ca and W sa , i = W ca , where T sa , i and W sa , i arethe temperature and humidity ratio of supply air to the i th zone. If a VAV box is not equipped with a reheat coil (cool-ing only), then the temperature of air supplied by it to its zonewill be at the conditioned air temperature, i.e., T sa , i = T ca .The control commands for a multi-zone VAV HVAC system with n z zones (i.e., VAV boxes) are: u : = ( m oa , T ca , m sa , i , T sa , i , i = , . . . , n z ) , (1)where m oa is the outdoor airflow rate, T ca is the conditionedair temperature, m sa , i is the supply airflow rate to the i th zone,and T sa , i is the supply air temperature to the i th zone. Notethat the humidity of conditioned air ( W ca ) which is suppliedto all the zones is indirectly controlled through T ca . Of the n z VAVs/zones in the building, n rhz VAVs are equipped witha reheat coil and n z − n rhz VAVs do not have a reheat coil(cooling only). For the latter, the supply air temperature willbe the same as the conditioned air temperature, i.e., T sa , i ( k ) = T ca ( k ) .The control commands in (1) are sent as set points tothe low-level control loops which are typically comprisedof proportional integral (PI) controllers. The role of a cli-mate control system is to vary these control commands sothat three main goals are satisfied: (i) ensure thermal com-fort, (ii) maintain indoor air quality, and (iii) use minimumamount of energy/cost.In an HVAC system as the one shown in Figure 1, thesupply duct pressure setpoint, p duct , is also usually a com-mand that the climate control system has to decide. We as-sume that the supply duct static pressure setpoint ( p duct ) iscontrolled based on “trim and respond” strategy [24], whichis commonly used in VAV systems, including in the Innova-tion Hub building that we use as a case study. The virtual building (VB) is a high-fidelity model of abuilding’s thermal dynamics and its HVAC system that willact as the plant for the controllers. The VB is chosen tomimic part of the Innovation Hub building in Gainesville,FL, USA, which is serviced by AHU-2 (among the twoAHUs that serve Phase I). Figures 2 and 3 show photos ofthe building and the relevant floor plans, respectively. Therooms supplied by the same VAV box are grouped together toform one large space (zone); the zones are enclosed by dot-ted lines in Figure 3. The first floor has 15 rooms which aregrouped into 9 zones, the second floor has 20 rooms whichare grouped into 12 zones, and the third floor has 21 roomswhich are grouped into 12 zones. In total, there are 56 roomsgrouped into 33 zones. The virtual building thus consistsof an air handling unit and 33 VAV boxes, of which 29 areequipped with reheat coils, and the remaining 4 do not havereheat coils (cooling only). The zones primarily consist ofoffices and labs.We use the Modelica library IDEAS (Integrated DistrictEnergy Assessment by Simulation) [25] to model the build-ing’s thermal dynamics.In order to model a zone we use the
RectangularZoneTemplate from the IDEAS li-brary. It consists of six components—which are a ceiling,a floor, and four walls—and an optional window. Thereare also external connections for each of the walls and theceiling. Depending on the usage, there are three types of3 hase-1 (a) Picture of the Innovation Hub building (view from south to north).Phase-1 is enclosed in the dashed lines.
Phase-1
AHU-2 c Google
NS EW (b) Top view of the Innovation Hub building. In this work, we consider airhandling unit 2 (AHU-2) which serves the southern half of Phase-1 (regionshaded in blue).
Fig. 2: Innovation Hub building located at the University ofFlorida campus.walls: (i) inner wall, which is used as a boundary betweenzones, (ii) outer wall, which is used as a boundary betweenoutside (atmosphere) and the zone, and (iii) boundary wall,which can be specified a fixed temperature or heat flow. Todefine a wall, dimensions, type of material, type of wall, andthe azimuth angle are required. The dimensions are obtainedfrom the mechanical drawings, the material type is chosenfrom the predefined materials available in the IDEAS library,the type of wall is chosen based on the zone’s location inthe building, and the azimuth angle is computed from thezone’s orientation. Windows are specified according to thedrawings, with the glazing material chosen from the IDEASlibrary. In this way, we model all the zones, which arethen connected appropriately to form the overall building;Figure 4 shows the model of floor 1. Since we are onlyinterested in modeling the southern half of Phase-1, thewalls that are adjacent to zones in the northern half areassumed to be at 22 . ◦ C (72 ◦ F ).Inputs to the building thermal dynamics portion of theVB are supply airflow rate ( m sa , i ), supply air temperature( T sa , i ), and supply air humidity ratio ( W sa , i ) for all the zones.These are implemented using the MassFlowSource T block from the IDEAS library; an ideal flow source thatproduces a specified mass flow with specified temperature,composition, and trace substances. Outputs of the simula-tor are temperature ( T z , i ) and humidity ratio ( W z , i ) of all thezones. The zone temperature and humidity are also influ-enced by several exogenous inputs: (i) outdoor weather con-ditions such as solar irradiation ( η sol ), outdoor air tempera-ture ( T oa ), etc. which are provided using the ReaderTMY3 block from the IDEAS library, (ii) internal sensible and la-tent heat loads due to occupants, which are computed based
Lunch & Vending room Conference Conference
Copy/Mailroom Men
Women Storage ITroom Mechanical
Storage
Telecom Electrical
Client room
IToffice
Recycling/Dumpsters
MechElevator
StairCorridor Corridor C o rr i d o r DoorWindow N VAV-104 VAV-105VAV-101 VAV-102 VAV-103 VAV-108
VAV-106 VAV-107
VAV-109 (a) Floor 1 plan.
Lab Lab LabLab LabLabOffice Office Office Office Office Office Office OfficeOffice Office
StorageTelecom ElectricalStorage Shared hood roomStorageElevator
StairDoorWindow C o rr i d o r C o rr i d o r Corridor Corridor N VAV-201 VAV-202 VAV-203 VAV-204VAV-206VAV-207 VAV-208 VAV-205SV-209 VAV-210 VAV-212VAV-211 (b) Floor 2 plan.
Office Office Office Office Office OfficeOffice Office Office OfficeLab Lab Lab Autoclave Lab Lab
Telecom
ElectricalElevator Storage
Shared Equipment
StairCorridor Corridor C o rr i d o r C o rr i d o r DoorWindow N VAV-301 VAV-304VAV-302 SV-303 VAV-306VAV-305 VAV-307VAV-312VAV-309 SV-310 VAV-311VAV-308 (c) Floor 3 plan.
Fig. 3: Floor plans of the southern half of Phase-1 which isserviced by AHU-2. v108
Flo B C D Cei A Flo B C D Cei A v107 Flo B C D Cei A Flo B C D Cei A v106 Flo B C D Cei A Flo B C D Cei A v105 Flo B C D Cei A Flo B C D Cei A v104 Flo B C D Cei A Flo B C D Cei A v102 Flo B C D Cei A Flo B C D Cei A v101 Flo B C D Cei A Flo B C D Cei A prescribedHeatFlow prescribedHeatFlow1prescribedHeatFlow2 prescribedHeatFlow3 prescribedHeatFlow4prescribedHeatFlow7prescribedHeatFlow8 sinWat i XiCm_ fl owT ms101XiCm_ fl owT ms108XiCm_ fl owT ms107XiCm_ fl owT ms106XiCm_ fl owT ms105XiCm_ fl owT ms104XiCm_ fl owT ms102XiCm_ fl owT ms109 v109 Flo B C D Cei A Flo B C D Cei A prescribedHeatFlow5x_pTphiTphiX_steam x_pTphi1TphiX_steamx_pTphi2TphiX_steam x_pTphi3TphiX_steamx_pTphi4TphiX_steam x_pTphi5TphiX_steamx_pTphi6TphiX_steamx_pTphi7TphiX_steamT_sa_108T_sa_107m_sa_108m_sa_107 W_z_108T_z_108T_z_107 W_z_107m_sa_106T_sa_106m_sa_105T_sa_105T_sa_104m_sa_104T_sa_102m_sa_102T_sa_101 T_z_101 W_z_101 W_z_102T_z_102T_z_104 W_z_104 W_z_105T_z_105 T_z_106 W_z_106q_other_101 co2_z_106co2_z_108co2_z_102co2_z_101 co2_z_105co2_z_104 co2_z_107q_other_102q_other_104 q_other_105 q_other_106q_other_108q_other_107n_101 n_108n_106n_105 n_107n_104 n_102W_sa ET_dryRHwind_speedT_dewT_skyDNIGHI m_sa_109T_sa_109 n_109T_z_109co2_z_109 W_z_109q_other_109co2 m_sa_101
12 34 5 SimInfoManager->ReaderTMY3 MassFlowSource_T RectangularZoneTemplate PrescribedHeatFlow X_pTphi Boundary_pT
Fig. 4: Floor 1 of the virtual building created in Dymola us-ing components from the IDEAS library [25].on the number of occupants provided to the zone block, and(iii) internal heat load due to lighting and equipment whichis given using the
PrescribedHeatFlow from the Mod-elica standard library.
Cooling and Dehumidifying Coil Model
The cooling coilmodel has five inputs and two outputs. The inputs are supplyairflow rate ( m sa ), mixed air temperature ( T ma ) and humidityratio ( W ma ), chilled water flow rate ( m w ), and chilled waterinlet temperature ( T wi ). The outputs are conditioned air tem-perature ( T ca ) and humidity ratio ( W ca ). We use a gray boxdata-driven model developed in our prior work [5]. The in-terested readers are referred to Section 2.1.2 of [5] for detailsregarding the model.4 ower Consumption Models For the HVAC system con-figuration presented in Figure 1, there are three main com-ponents which consume power. They are supply fan, coolingand dehumidifying coil, and reheating coils. The fan powerconsumption is modeled as: P f an ( k ) = α f an m sa ( k ) , (2)where m sa ( k ) is the total supply airflow rate at the AHU [26].The cooling and dehumidifying coil power consumptionis modeled to be proportional to the heat it extracts from themixed air stream: P cc ( k ) = m sa ( k ) (cid:2) h ma ( k ) − h ca ( k ) (cid:3) η cc COP c , (3)where h ma ( k ) and h ca ( k ) are the specific enthalpies of themixed and conditioned air respectively, η cc is the cooling coilefficiency, and COP c is the chiller coefficient of performance.Since a part of the return air is mixed with the outside air, thespecific enthalpy of the mixed air is: h ma ( k ) = r oa ( k ) h oa ( k ) + ( − r oa ( k )) h ra ( k ) , (4)where h oa ( k ) and h ra ( k ) are the specific enthalpies of the out-door and return air respectively, and r oa ( k ) is the outside airratio: r oa ( k ) : = m oa ( k ) m sa ( k ) . The specific enthalpy of moist airwith temperature T and humidity ratio W is given by [27]: h ( T , W ) = C pa T + W ( g H + C pw T ) , where g H is the heatof evaporation of water at 0 ◦ C , and C pa , C pw are specific heatof air and water at constant pressure.The reheating coil power consumption in the i th VAVbox is modeled to be proportional to the heat it adds to theconditioned air stream: P reheat , i ( k ) = m sa , i ( k ) C pa (cid:2) T sa , i ( k ) − T ca ( k ) (cid:3) η reheat COP h , (5)where η reheat is the reheating coil efficiency, and COP h is theboiler coefficient of performance. Overall Plant
The overall plant, i.e., virtual building—consisting of the building thermal model, cooling and dehu-midifying coil model, and power consumption models—issimulated using SIMULINK and MATLAB©. The build-ing thermal model is constructed in DYMOLA 2021 andit is exported into an FMU (Functional Mockup Unit). Itis then imported into SIMULINK using the FMI Kit forSIMULINK. The remaining models are constructed directlyin SIMULINK.
MZHC ) Recall that both the proposed and the baseline con-trollers need to decide the following control commands: u ( k ) : = [ m oa ( k ) , T ca ( k ) , m sa , i ( k ) , T sa , i ( k )] T ∈ ℜ + n z + n rhz . Figure 5 shows the structure of the proposed
MZHC . Thehigh-level controller is based on MPC and decides the con-trol commands for the AHU: outdoor air flow rate ( m oa ) andconditioned air temperature ( T ca ). The low-level controller isa projection-based feedback controller and decides the con-trol commands for each of the VAV boxes/zones: supply airflow rate ( m sa , i ) and supply air temperature ( T sa , i ). Thesecontrollers are described in detail next. The high-level controller (HLC) is based on MPC thatuses a low-resolution model of the building which is dividedinto a small number of meta-zones. Each meta-zone is an ag-gregation of multiple zones in the real building. This aggre-gation can be done in any number of ways. In this work weaggregate all the zones in a floor into a meta-zone, which isdenoted by f ∈ F : = { , . . . , n f } , where n f is the total numberof floors/meta-zones. The Innovation Hub building has threefloors, so we aggregate it into three meta-zones. The set ofall VAVs/zones in floor f is denoted as I f (so | ∪ f ∈ F I f | = n z ),of which those equipped with reheat coils is denoted as I rh , f (so | ∪ f ∈ F I rh , f | = n rhz ). The HLC decides on the followingcontrol commands based on the aggregate models: u HLC ( k ) : = (cid:0) m oa ( k ) , T ca ( k ) , m sa , f ( k ) , T sa , f ( k ) , ∀ f ∈ F (cid:1) ∈ ℜ +( × n f ) , (6)where m sa , f ( k ) : = ∑ i ∈ I f m sa , i ( k ) is the aggregate (total) supplyairflow rate to all the zones in floor/meta-zone f and T sa , f ( k ) is the aggregate supply air temperature. Of the control com-mands computed in (6), m oa ( k ) and T ca ( k ) can be directlysent to the plant. The remaining information computed bythe HLC including m sa , f ( k ) and T sa , f ( k ) are used by the low-level controller (LLC), described in Section 3.2, to decide on Tcamoa Tz,f
Building+HVAC SystemHigh-LevelMPCLow-LevelProjectionControllerState Estimator
Tsa,imsa,imsaPreheat W z,f
Tz,iTz,f Tw,f[qint,f, ω int,f] Tz,fHLC WeatherForecast
HLCHLC
OutputsControl Commands
TcaHLC
Fig. 5: Structure of the proposed multi-zone hierarchicalcontroller (
MZHC ). We denote estimates as ˆ • and forecastsas ˆˆ • . Variables with a subscript i are for the individual zones,while the variables with a subscript f represent the aggregatequantities for each floor/meta-zone. In this figure, T z , f , W z , f , q ac , f , ˆ T z , f , ˆ T w , f , ˆˆ q int , f , ˆˆ ω int , f , and T HLCz , f are ∀ f ∈ F ; m sa , i and T z , i are for i ∈ I f , ∀ f ∈ F ; T sa , i is for i ∈ I rh , f , ∀ f ∈ F .5he supply airflow rate ( m sa , i ( k ) ) and supply air temperature( T sa , 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 asubscript f represent the aggregate quantities for each meta-zone.For MPC formulation, we use a model interval of ∆ t = ∆ T =
15 minutes, and aprediction/planning interval of T =
24 hours. So we have T = N ∆ T and ∆ T = M ∆ t , where N =
96 (planning hori-zon) and M =
3. The control inputs for N time steps areobtained by solving an optimization problem of minimizingthe energy consumption subject to thermal comfort, indoorair quality, and actuator constraints. Then the control com-mands obtained for the first time step are sent to the plantand the LLC. The optimization problem is solved again forthe next N time steps with the initial states of the model ob-tained from a state estimator, which uses measurements fromthe plant. This process is repeated at the next control timestep, i.e., after an interval of ∆ T .To describe the optimization problem, first we define thestate vector x ( k ) and the vector of control commands andinternal variables v ( k ) as: x ( k ) : = (cid:0) T z , f ( k ) , T w , f ( k ) , W z , f ( k ) , ∀ f ∈ F (cid:1) ∈ ℜ × n f , (7) v ( k ) : = (cid:0) u HLC ( k ) , m w ( k ) , W ca ( k ) (cid:1) ∈ ℜ +( × n f )+ , (8)where T z , f ( k ) , T w , f ( k ) , and W z , f ( k ) are the aggregatezone temperature, wall temperature, and humidity ratio offloor/meta-zone f , respectively; u HLC ( k ) is the control com-mand vector defined in (6) and m w ( k ) is the chilled waterflow rate into the cooling coil. The exogenous input vectoris: w ( k ) : = (cid:0) η sol ( k ) , T oa ( k ) , W oa ( k ) , q int , f ( k ) , ω int , f ( k ) , ∀ f ∈ F (cid:1) ∈ ℜ +( × n f ) , (9)where η sol ( k ) is the solar irradiance, T oa ( k ) is the out-door air temperature, W oa ( k ) is the outdoor air humid-ity ratio, q int , f ( k ) is the aggregate internal heat load infloor/meta-zone f due to occupants, lights, equipment, etc.,and ω int , f ( k ) is the aggregate rate of water vapor gen-eration in floor/meta-zone f due to occupants and othersources. We denote the forecast of these exogenous inputsas ˆˆ w ; in Section 5, we discuss how these forecasts are ob-tained. The vector of nonnegative slack variables ζ ( k ) : = (cid:0) ζ lowT , f ( k ) , ζ highT , f ( k ) , ζ lowW , f ( k ) , ζ highW , f ( k ) , ∀ f ∈ F (cid:1) ∈ ℜ × n f , is in-troduced for feasibility of the optimization problem.The optimization problem at time index j is:min V , X , Z j + NM − ∑ k = j (cid:20) P f an ( k ) + P cc ( k ) + ∑ f ∈ F P reheat , f ( k ) (cid:21) ∆ t + P slack ( k ) , (10a) where P f an ( k ) is given by (2), P cc ( k ) is givenby (3), P reheat , f ( k ) : = m sa , f ( k ) C pa (cid:2) T sa , f ( k ) − T ca ( k ) (cid:3) η reheat COP h , V : = [ v T ( j ) , v T ( j + ) , . . . , v T ( j + NM − )] T , X : = [ x T ( j + ) , x T ( j + ) , . . . , x T ( j + NM )] T , and Z : = [ ζ T ( j + ) , ζ T ( j + ) , . . . ζ T ( j + NM )] T . The lastterm, P slack , penalizes the aggregate zone temperature andhumidity slack variables: P slack ( k ) : = ∑ f ∈ F (cid:20) λ lowT ζ lowT , f ( k + ) + λ highT ζ highT , f ( k + )+ λ lowW ζ lowW , f ( k + ) + λ highW ζ highW , f ( k + ) (cid:21) , where the λ s are penalty parameters. The total supply air-flow rate m sa ( k ) used in P f an ( k ) and P cc ( k ) , is given by m sa ( k ) = ∑ f ∈ F m sa , f ( k ) = ∑ f ∈ F ∑ i ∈ I f m sa , i ( k ) . The optimal controlcommands are obtained by solving the optimization prob-lem (10a) subject to the following constraints: T z , f ( k + ) = T z , f ( k ) + ∆ t (cid:20) T oa ( k ) − T z , f ( k ) τ za , f + T w , f ( k ) − T z , f ( k ) τ zw , f + A z , f η sol ( k ) + q int , f ( k ) + q ac , f ( k ) C z , f (cid:21) (10b) T w , f ( k + ) = T w , f ( k ) + ∆ t (cid:20) T oa ( k ) − T w , f ( k ) τ wa , f + T z , f ( k ) − T w , f ( k ) τ wz , f + A w , f η sol ( k ) (cid:21) (10c) W z , f ( k + ) = W z , f ( k ) + ∆ tR g T z , f ( k ) V f P da (cid:20) ω int , f ( k )+ m sa , f ( k ) W ca ( k ) − W z , f ( k ) + W ca ( k ) (cid:21) (10d) T ca ( k ) = T ma ( k ) + m w ( k ) f (cid:0) T ma ( k ) , W ma ( k ) , m sa ( k ) , m w ( k ) (cid:1) (10e) W ca ( k ) = W ma ( k ) + m w ( k ) g (cid:0) T ma ( k ) , W ma ( k ) , m sa ( k ) , m w ( k ) (cid:1) (10f) T lowz , f ( k ) − ζ lowT , f ( k ) ≤ T z , f ( k ) ≤ T highz , f ( k ) + ζ highT , f ( k ) (10g) a low T z , f ( k ) + b low − ζ lowW , f ( k ) ≤ W z , f ( k ) ≤ a high T z , f ( k ) + b high + ζ highW , f ( k ) (10h) m minoa ≤ m oa ( k ) ≤ m maxoa (10i) T ca ( k + ) ≤ min (cid:0) T ca ( k ) + T rateca ∆ t , T ma ( k + ) , T highca (cid:1) (10j) T ca ( k + ) ≥ max (cid:0) T ca ( k ) − T rateca ∆ t , T lowca (cid:1) (10k) W ca ( k ) ≤ W ma ( k ) (10l) r oa ( k ) = m oa ( k ) / m sa ( k ) (10m) r oa ( k + ) ≤ min (cid:0) r oa ( k ) + r rateoa ∆ t , r highoa (cid:1) (10n) r oa ( k + ) ≥ max (cid:0) r oa ( k ) − r rateoa ∆ t , r lowoa (cid:1) (10o)6 lowsa , f ≤ m sa , f ( k ) ≤ m highsa , f (10p) T ca ( k ) ≤ T sa , f ( k ) ≤ T highsa (10q) ζ lowT , f ( k + ) , ζ highT , f ( k + ) , ζ lowW , f ( k + ) , ζ highW , f ( k + ) ≥ ∀ f ∈ F . Constraints (10b)-(10f), (10i), (10l)-(10m), and(10p)-(10r) are for k ∈ { j , . . . , j + NM − } , constraints (10g)and (10h) are for k ∈ { j + , . . . , j + NM } , and constraints(10j)-(10k) and (10n)-(10o) are for k ∈ { j − , . . . , j + NM − } . The control commands remain the same for M time steps,as the control interval ∆ T = M ∆ t , i.e., u HLC ( k ) = u HLC ( k + ) = · · · = u HLC ( k + M − ) , ∀ k ∈ { j , j + M , . . . , j + NM − } .Constraints (10b) and (10c) are due to the aggregatethermal dynamics of floor/meta-zone f , which is a dis-cretized form of an RC (resistor-capacitor) network model,specifically a 2R2C model. The two states of the modelare aggregate zone temperature ( T z , f ) and wall temperature( T w , f , a fictitious state). In constraint (10b), q ac , f ( k ) isthe heat influx due to the HVAC system and is given by q ac , f ( k ) : = m sa , f ( k ) C pa ( T sa , f ( k ) − T z , f ( k )) . The model hasseven parameters { C z , f , τ zw , f , τ za , f , A z , f , τ wz , f , τ wa , f , A w , f } . Inthe evaluation study, they are estimated using the algorithmpresented in [20] and will be discussed later in Section 5.The constraint (10d) is for the aggregate humidity dy-namics of floor/meta-zone f , where W z , f is the aggregatezone humidity ratio, V f is the volume of meta-zone f , R g isthe specific gas constant of dry air, P da is the partial pressureof dry air, and W ca is the conditioned air humidity ratio [28].Constraints (10e) and (10f) are for the control-orientedcooling and dehumidifying coil model, which was developedin our prior work [5]. The specific functional form in (10e)and (10f) is chosen so that when the chilled water flow rate iszero, no cooling or dehumidification of the air occurs, so theconditioned air temperature and humidity ratio are equal tothe mixed air temperature and humidity ratio: T ca = T ma and W ca = W ma , when m w =
0. The interested readers are referredto [5, Section 3.1.1] for details regarding the model.Constraints (10g) and (10h) are box constraints to main-tain the temperature and humidity of the meta-zones withinthe allowed comfort limits. The constraints are softened us-ing slack variables ζ lowT , f ( k ) , ζ highT , f ( k ) , ζ lowW , f ( k ) , and ζ highW , f ( k ) ;constraint (10r) ensures that these slack variables are non-negative. Imposing constraints directly on the relative hu-midity of zones ( RH z ) is difficult, as relative humidity is ahighly nonlinear function of dry bulb temperature and hu-midity ratio [27, Chapter 1]. So we linearize this functionwhich gives us a low , b low , a high , and b high in (10h), and helpsin converting the constraints on relative humidity to humidityratio ( W z ).Constraint (10i) is for the outdoor airflow rate, wherethe minimum allowed value ( m minoa ) is computed based on theventilation requirements specified in ASHRAE 62.1 [29] andto maintain positive building pressurization.Constraints (10j)-(10k) and (10n)-(10o) are to take intoaccount the capabilities of the cooling coil and damper actu- ators. In constraints (10j) and (10l), the inequalities T ca ( k + ) ≤ T ma ( k + ) and W ca ( k ) ≤ W ma ( k ) ensure that the coolingcoil can only cool and dehumidify the mixed air stream; itcannot add heat or moisture. Similarly, in constraint (10q)the inequality T sa , f ( k ) ≥ T ca ( k ) ensures that the reheat coilscan only add heat; it cannot cool.Constraint (10p) is to take into account the capabilitiesof the fan and aggregate capabilities of the VAV boxes. Thelimits m lowsa , f and m highsa , f are computed using the VAV sched-ule from the mechanical drawings of a building as follows: m lowsa , f : = ∑ i ∈ I f m lowsa , i and m highsa , f : = ∑ i ∈ I f m highsa , i .Note that of the states x ( k ) defined in (7), T w , f is a ficti-tious state that cannot be measured, while the other two statesaggregate zone temperature ( T z , f ) and aggregate zone humid-ity ratio ( W z , f ) are measured. So we estimate the current stateˆ x ( k ) using a Kalman filter. The role of the low-level controller (LLC) is to ap-propriately distribute the aggregate quantities—such as thetotal supply airflow rate and reheat power consumption—computed by the HLC to individual zones/VAVs. The LLCneeds to do so by capturing two important properties: (i) itshould consider the needs of individual zones and distributeaccordingly, and (ii) it should act in coherence with the HLC,so that there is minimal mismatch for the MPC optimizationin the next round.The LLC is a projection-based feedback controller thatdecides on the supply airflow rate and supply air temperaturefor each VAV box/zone. That is, the control command vectorthat the LLC needs to decide is: u LLC ( k ) : = [ m sa , i ( k ) , T sa , i ( k )] T ∈ ℜ n z + n rhz , where for m sa , i , i ∈ I f , ∀ f ∈ F , and for T sa , i , i ∈ I rh , f , ∀ f ∈ F .It decides these control commands by using the following in-formation from the HLC: (i) total allowed supply airflow rateto all the zones m sa ( k ) = ∑ f ∈ F m sa , f ( k ) , (ii) total allowed re-heat power consumption P reheat ( k ) = ∑ f ∈ F P reheat , f ( k ) , (iii) thetemperature at which the zones in each meta-zone should bemaintained at T z , f ( k + ) , and (iv) the conditioned air tem-perature T ca ( k ) . Here on in this section, we will be using thesuperscript HLC ( • HLC ) for these variables to make it clearthat these are obtained from the high-level controller.First the needs of each zone are assessed based onthe current measured temperature T z , i ( k ) and the range itshould be in [ T htgz , i ( k ) , T clgz , i ( k )] and are translated into the de-sired supply airflow rate m dsa , i ( k ) and supply air temperature T dsa , i ( k ) . Then these desired values along with the informa-tion obtained from the HLC are used to solve a projectionproblem to compute the control commands for all the zones, u LLC ( k ) .The procedure used to compute the desired values m dsa , i ( k ) and T dsa , i ( k ) is explained below. This is similar to7he Dual Maximum control logic presented in Section 4; aschematic representation of it is shown in Figure 6.1. First the temperature range [ T htgz , i ( k ) , T clgz , i ( k )] in whicheach zone should be is computed as follows: T htgz , i ( k ) = max (cid:0) T HLCz , f ( k + ) − ˜ T dbz / , T lowz , f (cid:1) ∀ i ∈ I f and T clgz , i ( k ) = min (cid:0) T HLCz , f ( k + ) + ˜ T dbz / , T highz , f (cid:1) ∀ i ∈ I f ,where T HLCz , f ( k + ) is obtained from the HLC, ˜ T dbz isa deadband, and T lowz , f and T highz , f are the limits used inconstraint (10g).2. If the zone temperature is between the cooling and heat-ing setpoints ( T z , i ( k ) ∈ [ T htgz , i ( k ) , T clgz , i ( k )] ), then the con-troller is in deadband mode. The supply airflow rate isdesired to be at its minimum and no heating is required,i.e., m dsa , i ( k ) = m lowsa , i and T dsa , i ( k ) = T HLCca ( k ) .3. If the zone temperature is warmer than the cooling set-point ( T z , i ( k ) > T clgz , i ( k ) ), then the controller is in coolingmode. The supply airflow rate is desired to be increasedas needed and no heating is required, i.e., m dsa , i ( k ) = min (cid:0) m lowsa , i + K clgm , i ( T z , i ( k ) − T clgz , i ( k )) , m highsa , i (cid:1) and T dsa , i ( k ) = T HLCca ( k ) .4. If the zone temperature is cooler than the heating set-point ( T z , i ( k ) < T htgz , i ( k ) ), then the controller is in heatingmode. Heating is required, and the supply airflow rateis desired to be increased only if additional heating isneeded, i.e., T dsa , i ( k ) = min (cid:0) T HLCca ( k ) + K htgT , i ( T htgz , i ( k ) − T z , i ( k )) , T highsa (cid:1) ; if T dsa , i ( k ) = T highsa , then m dsa , i ( k ) = min (cid:0) m lowsa , i + K htgm , i ( T htgz , i ( k ) − T z , i ( k )) , m high , reheatsa , i (cid:1) , other-wise m dsa , i ( k ) = m lowsa , i .5. Finally, we impose the following rate constraints: m sa , i ( k − M ) − m ratesa , i ∆ T ≤ m dsa , i ( k ) ≤ m sa , i ( k − M ) + m ratesa , i ∆ T , T sa , i ( k − M ) − T ratesa , i ∆ T ≤ T dsa , i ( k ) ≤ T sa , i ( k − M ) + T ratesa , i ∆ T , where m sa , i ( k − M ) and T sa , i ( k − M ) are the supply air-flow rate and supply air temperature from the previouscontrol time step.These desired values— m dsa , i ( k ) and T dsa , i ( k ) —along withinformation from the HLC are used to solve the followingprojection problem to obtain the control commands for allthe zones, u LLC ( k ) :min u LLC ( k ) ∑ f ∈ F ∑ i ∈ I f λ m , i ( m sa , i ( k ) − m dsa , i ( k )) + ∑ f ∈ F ∑ i ∈ I rh , f λ T , i ( T sa , i ( k ) − T dsa , i ( k )) (11a) subject to the following constraints: ∑ f ∈ F ∑ i ∈ I f m sa , i ( k ) ≤ m HLCsa ( k ) (11b) ∑ f ∈ F ∑ i ∈ I rh , f m sa , i ( k ) C pa (cid:16) T sa , i ( k ) − T HLCca ( k ) (cid:17) η reheat COP h ≤ P HLCreheat ( k ) (11c) m lowsa , i ≤ m sa , i ( k ) ≤ m highsa , i , ∀ i ∈ I f , ∀ f ∈ F (11d) T ca ( k ) ≤ T sa , i ( k ) ≤ T highsa , ∀ i ∈ I rh , f , ∀ f ∈ F (11e)where the sets I f and I rh , f are defined at the beginning ofthis section, λ s are weights, m HLCsa ( k ) = ∑ f ∈ F m HLCsa , f ( k ) , and P HLCreheat ( k ) = ∑ f ∈ F P HLCreheat , f ( k ) .Constraints (11b) and (11c) are to ensure that the totalsupply airflow rate and reheat power consumption do not ex-ceed the limits computed by the HLC. Constraints (11d) and(11e) are to take in to account the capabilities of the VAVboxes and reheat coils. In constraint (11e), the inequality T sa , i ( k ) ≥ T ca ( k ) ensures that the reheat coils can only addheat to the conditioned air and cannot cool. The upper limiton supply air temperature ( T highsa ) in constraint (11e) is to pre-vent thermal stratification [22]. BL ) For zone climate control, we consider the
Dual Maxi-mum algorithm [22] as the baseline; a schematic representa-tion of this algorithm is shown in Figure 6. Even though
Sin-gle Maximum is more commonly used, including in the Inno-vation Hub building, we choose
Dual Maximum for the base-line, as it is more energy-efficient among the two [22, 30].The
Dual Maximum controller operates in three modes basedon the zone temperature ( T z , i ): (i) cooling, (ii) heating, and(iii) deadband. The zone’s supply airflow rate ( m sa , i ) andsupply air temperature ( T sa , i ) are varied based on the mode,as explained below.1. Cooling mode:
If the zone temperature is warmer thanthe cooling setpoint, then the controller is in coolingmode. The supply airflow rate ( m sa , i ) is varied betweenthe minimum ( m lowsa , i ) and maximum ( m highsa , i ) as needed, Maximum supply airtemperature
Supply airtemperature (T sa,i ) Maximum cooling airflow rateSupply airflow rate (m sa,i )Deadband mode Cooling modeHeating mode
ZonetemperatureSupply airflow rate
Minimumairflow rate
Supply airtemperature
Heatingsetpoint Cooling setpoint
Maximum heatingairflow rate
Fig. 6: Schematic of the
Dual Maximum control algorithm.8nd the supply air temperature ( T sa , i ) is equal to the con-ditioned air temperature ( T ca ), i.e., no reheat.2. Heating mode:
If the zone temperature is below theheating setpoint, then the controller is in heating mode.First, the supply air temperature ( T sa , i ) is increased upto the maximum allowed value ( T highsa ) as needed tomaintain the zone temperature at the heating setpoint.If the zone temperature still cannot be maintained atthe heating setpoint, then the supply airflow rate is in-creased between the minimum ( m lowsa , i ) and the heatingmaximum ( m high , reheatsa , i ) values.3. Deadband mode:
If the zone temperature is betweenthe heating and cooling setpoints, then the controller isin deadband mode. The supply airflow rate is kept at theminimum, and the supply air temperature is equal to theconditioned air temperature, i.e., no reheat.In the case of VAV boxes that do not have reheat coils, thelogic during cooling and deadband modes are the same. Inheating mode, the supply airflow rate is at the minimum andthe supply air temperature is equal to the conditioned air tem-perature, as the VAV box cannot heat.For the AHU-level commands, the BL controller usesfixed conditioned air temperature that is determined based onexpected thermal (sensible and latent) load, and fixed out-door airflow rates based on ventilation requirements, e.g.,ASHRAE 62.1 [29]. Another consideration in choosing out-door airflow rate is building positive pressurization require-ments [27]. Recall that the plant is based on an air handling unitserving 33 zones, of which 29 are equipped with reheatcoils, and the remaining 4 do not have reheat coils (cool-ing only). See Table 1 and Figure 3 for the entire list of VAVboxes/zones. Of the 29 VAV boxes with reheat, three of themserve laboratories which are equipped with fume hoods (209,303, and 310), and one of them serve restrooms (103). TheVAV boxes serving these labs need to be controlled to sat-isfy the negative pressurization requirements with respect tocorridor, so we assume that they operate according to theexisting rule based feedback control strategy. Therefore, n z =
29 and n rhz =
25; for m sa , i , i ∈ { I : = { , } , I : = { , } , I : = { , } and for T sa , i , i ∈ { I rh , : = { I \{ , }} , I rh , : = { I \{ }} , I rh , : = { I \{ }}} . The sets I , I , and I defined above are the VAVs/zones in floors 1, 2, and 3, re-spectively. The sets I rh , , I rh , , and I rh , exclude the VAVboxes which do not have a reheat coil.The outdoor weather data used in simulations is ob-tained from the National Solar Radiation Database [31] forGainesville, Florida. As mentioned in Section 2.1, the in-ternal heat load due to occupants are computed based on thenumber of occupants provided to the zone block. We assumethat the zones are occupied from Monday to Friday between8:00 a.m. to noon and 1:00 p.m. to 5:00 p.m., with the to-tal number of occupants ( n p , f ) in floor 1 as 24, in floor 2 Table 1: VAV schedule. VAV Reheat m lowsa , i ( kg / s ) m high , reheatsa , i ( kg / s ) m highsa , i ( kg / s )
101 Yes 0.27 0.68 1.36102 Yes 0.05 0.10 0.20103 Yes 0.07 0.17 0.34104 Yes 0.57 0.57 1.14105 Yes 0.30 0.30 0.60106 No 0.04 - 0.23107 No 0.13 - 0.45108 Yes 0.14 0.33 0.66109 Yes 0.08 0.21 0.39201 Yes 0.21 0.52 1.03202 Yes 0.11 0.28 0.57203 Yes 0.11 0.28 0.57204 Yes 0.06 0.14 0.28205 No 0.13 - 0.57206 Yes 0.07 0.16 0.33207 Yes 0.13 0.32 0.64208 Yes 0.09 0.20 0.41209 Yes 0.18 0.18 0.57210 Yes 0.11 0.28 0.57211 Yes 0.10 0.26 0.51212 Yes 0.14 0.34 0.68301 Yes 0.16 0.41 0.82302 Yes 0.11 0.28 0.57303 Yes 0.28 0.28 0.57304 Yes 0.13 0.32 0.64305 Yes 0.10 0.23 0.47306 Yes 0.11 0.28 0.57307 No 0.13 - 0.57308 Yes 0.11 0.28 0.57309 Yes 0.11 0.28 0.57310 Yes 0.28 0.28 0.57311 Yes 0.16 0.40 0.80312 Yes 0.10 0.26 0.51 as 26, and in floor 3 as 22. We assume a power density of12.92 W/m (1.2 W/ft ) for internal heat load due to light-ing and equipment. For special purpose rooms like electricaland telecommunication, we use a higher power density of53.82 W/m (5 W/ft ). These heat loads from lighting andequipment are assumed to be halved during weekends.The following zone temperature and humidity limitsare used in the simulations: T lowz = 21 . ◦ C (70 ◦ F), T highz =23 . ◦ C (74 ◦ F), RH lowz = 20%, and RH highz = 65%. The cho-9 omfort envelope Fig. 7: Thermal comfort envelope from [27] shown as thehatched areas. Comfort envelope chosen in this paper shownas the green shaded area.sen thermal comfort envelope is shown in Figure 7. Typicallythe zone temperature limits during unoccupied mode (unocc)are relaxed when compared to the occupied mode (occ), i.e., [ T low , occz , T high , occz ] ⊆ [ T low , unoccz , T high , unoccz ] . Due to its us-age, the Innovation Hub building is always operated in oc-cupied mode, so we assume the same in simulations. Forthe simulation results reported later, the zone temperature vi-olation is computed as max (cid:0) T z ( k ) − T highz , T lowz − T z ( k ) , (cid:1) and the zone relative humidity violation is computed as max (cid:0) RH z ( k ) − RH highz , RH lowz − RH z ( k ) , (cid:1) , with the upperand lower limits mentioned above.The fan power coefficient α f an in (2) is14.2005 W / ( kg / s ) , which is obtained using a leastsquares fit to data collected from the building. The parame-ters of the cooling and dehumidifying coil model used in theplant are fit using the procedure explained in Section 2.1.2of [5]. The root mean square errors for the validation dataset are 0.25 ◦ C (0.46 ◦ F, 2%) for T ca and 0.22 × − kg w / kg da (2.6%) for W ca .The AHU in the building is equipped with a draw-through supply fan and therefore the fan is located after thecooling coil. The fan emits heat, which leads to a slight in-crease in the conditioned air temperature before it is suppliedto the VAV boxes. For the simulations, we assume this in-crease in temperature to be 1.11 ◦ C (2 ◦ F).
MZHC Parameters:
The optimization problemsin HLC and LLC are solved using CasADi [32] andIPOPT [33], a nonlinear programming (NLP) solver, on aDesktop Windows computer with 16GB RAM and a 3.60GHZ × ∆ t = ∆ T =
15 minutes, T =
24 hours, N =
96, and M = Parameter Units Floor 1 Floor 2 Floor 3 C z , f kWh / o C τ zw , f hours 0.5108 2.2161 2.5622 τ za , f hours 200 150 150 A z , f o Cm / kWh τ wz , f hours 18.7779 68.4388 100 τ wa , f hours 4157.5 1145.7 1129.2 A w , f o Cm / kWh × − × − × − Time (Hours)
Fig. 8: Out of sample aggregate zone temperature ( T z , f ) pre-diction results using the estimated aggregate RC networkmodel.are 0.97 ◦ C (1.75 ◦ F, 7.6%) for T ca and 0.63 × − kg w / kg da (7.6%) for W ca .Since the Innovation Hub building has three floors, weaggregate it into three meta-zones, i.e., f ∈ { , , } = : F .The parameters of the aggregate thermal dynamics model foreach meta-zone are estimated using the algorithm presentedin [20]. The parameters are shown in Table 2. Figure 8 showsthe out of sample prediction results using the estimated ag-gregate RC network model.For the aggregate humidity dynamics model, floor vol-umes used are V = . m , V = . m , and V = . m , which are obtained from mechanical drawings ofthe building.The following limits are used for the zone tem-perature constraint (10g): T lowz , f = ◦ C (70 ◦ F) and T highz , f = ◦ C (74 ◦ F). The coefficients for the humidityconstraint in (10h) are a high = . kg w / kg da / ◦ C and b high = − . kg w / kg da , which corresponds to a rela-tive humidity of 60%, and a low = . kg w / kg da / ◦ C and b low = − . kg w / kg da , which corresponds to a relativehumidity of 20%. We introduce a factor of safety by using aslightly tighter higher limit of 60% for the relative humidityof the zones when compared to the thermal comfort envelope10resented in Figure 7.The minimum allowed value for the outdoor airflowrate ( m minoa ) is 3.24 kg/s (5700 cfm), which is obtained fromthe AHU schedule in the mechanical drawings for the build-ing. The maximum possible value for the outdoor airflowrate ( m maxoa ) is 8.52 kg/s (15000 cfm). The various lim-its on the supply airflow rates are obtained using the VAVschedule presented in Table 1. The remaining limits usedin the controllers are as follows: r lowoa = r highoa = T dbz = . ◦ C (1 ◦ F), T lowca = . ◦ C (53 ◦ F), T highca = . ◦ C(63 ◦ F), and T highsa = ◦ C (86 ◦ F). The higher limit on theconditioned air temperature ( T highca ) is to introduce a factor ofsafety and make the controller robust.The MPC controller requires predictions of the variousexogenous inputs specified in (9). We compute the loads dueto occupants in q int , f and ω int , f based on the occupancy pro-file used in simulating the plant. The outdoor weather relatedexogenous inputs are assumed to be fully known. BL Parameters:
The cooling and heating setpoints arechosen to be 21.1 ◦ C (70 ◦ F) and 23.3 ◦ C (74 ◦ F), respectively.The minimum, maximum, and heating maximum values forthe supply airflow rate of all the VAV boxes are listed in Ta-ble 1. The maximum allowed value for the supply air tem-perature ( T highsa ) is 30 ◦ C (86 ◦ F). The conditioned air temper-ature ( T ca ) is kept at a constant value of 11.67 ◦ C (53 ◦ F). Typ-ically T ca is kept at 12.8 ◦ C (55 ◦ F), especially in hot-humidclimates, to ensure humidity control [34], but recall that weassume there is a 1.11 ◦ C (2 ◦ F) increase in temperature be-cause of the heat from the draw-through supply fan in theAHU, so we keep it at 11.67 ◦ C (53 ◦ F) to compensate forit. The outdoor airflow rate is kept at 3.24 kg/s (5700 cfm),which is obtained from mechanical drawings for the build-ing.
Performance of the controllers is compared using threetypes of outdoor weather conditions: hot-humid (Jul/06 toJul/13), mild (Feb/19 to Feb/26), and cold (Jan/30 to Feb/06).
BL Hot-HumidMZHC BL MildMZHC BL ColdMZHC E t o t a l f o r a w ee k ( k W h )
89% 38% 32%
FanCoolingReheat
Fig. 9: Comparison of the total energy consumed for a weekwhen using the baseline ( BL ) and proposed ( MZHC ) con-trollers for different outdoor weather conditions. The proposed controller reduces energy use significantlycompared to BL , from approximately 11% to 68% dependingon weather; see Figure 9. The indoor climate control perfor-mance of MZHC and BL are nearly identical. With MZHC ,the RMSE (root mean square error) of zone temperature vi-olation is 0.1 ◦ C (0.18 ◦ F) and the RMSE of zone relative hu-midity violation is 0.05%, while with BL they are 0.01 ◦ C(0.02 ◦ F) and 0% respectively. The computational cost of theproposed
MZHC is small, just a few seconds are needed tocompute decisions at every control update. On an average ittakes 3.28 seconds to solve the optimization problem in theHLC and 0.018 seconds to solve the optimization problem inthe LLC.Simulation results for the different weather conditionsare discussed in detail next.
Figure 10 shows the simulation time traces for a hot-humid week. It is found that using the proposed
MZHC leadsto 11% energy savings when compared to BL , as presentedin Figure 9.Both controllers lead to negligible violations of aggre-gate zone temperature ( T z , f ) and relative humidity ( RH z , f )constraints. Data for the three meta-zones are shown in Fig-ure 10(c), and for three individual zones, one from each floor,are shown in Figure 11. BL ensures that dry enough air issupplied to the zones at all times by keeping the conditionedair temperature ( T ca ) at a constant low value of 11.67 ◦ C(53 ◦ F), and hence the humidity limit is not violated. In thecase of
MZHC , the humidity constraint is found to be ac-tive always. This can be seen in Figure 10(c); recall that weuse a tighter constraint of 60% instead of 65% to introduce afactor of safety. This active constraint limits the increase in T ca , which can be seen in Figure 10(d). One of the reasonsreheating is required even during such a hot week (the out-door air temperature is as high as 32.2 ◦ C/90 ◦ F) is because ofthis active humidity constraint, which requires dry, and thus,cold air to be supplied to the zones. This could also be oneof the reasons
MZHC decides to maintain the zone tempera-tures at the lower limit (see Figure 10(c)). Keeping the zonesat a higher temperature will require an increase in reheatingenergy.Most of the prior works on using MPC for HVAC con-trol ignore humidity and latent heat in their formulation. Inan attempt to reduce energy/cost, such controllers are likelyto make decisions during these hot-humid conditions whichwill lead to poor indoor humidity, as they are unaware of thefactors mentioned above [5]. Thus such controllers cannotbe used particularly in hot-humid climates.The energy savings by
MZHC is because of two mainreasons. One,
MZHC increases the T ca as long as the hu-midity constraints are not violated, while BL uses a conser-vatively designed value as explained above. This leads toa reduction in cooling energy consumption by MZHC ; seeFigures 9 and 10(b). Two, the warmer T ca supplied to theVAV boxes requires lesser reheating in the case of MZHC .This leads to a reduction in the reheat energy consumption,11
24 48 72 96 120 144 168202530350 24 48 72 96 120 144 1680501000 24 48 72 96 120 144 168
Time (Hours) (a) Outdoor weather data used in simulations (outdoor air temperature, out-door air relative humidity, and solar irradiance).
MZHCBL
24 48 72 96 120 144 168025500 24 48 72 96 120 144 168
Time (Hours) (b) Comparing the power consumptions (fan, cooling, and total reheatpower).(c) Aggregate conditions of floors/meta-zones 1, 2, and 3 (temperatures andrelative humidities) when using
MZHC and BL . The black dashed lines arethe thermal comfort limits. MZHCBL
Time (Hours) (d) Conditions at the AHU (supply airflow rate, outdoor airflow rate, con-ditioned air temperature, and conditioned air humidity ratio).
Fig. 10: Comparison of the two controllers for a hot-humid week (Jul/06 to Jul/13, Gainesville, Florida, USA).Fig. 11: Individual zone conditions (temperatures and rela-tive humidities) when using
MZHC and BL for a hot-humidweek. The black dashed lines are the thermal comfort limits. which can be seen in Figures 9 and 10(b). The decisions re-garding the supply airflow rates are found to be the same forboth the controllers, and thus the fan energy consumptionsare identical.Since the outdoor air is hot and humid (see Fig-ure 10(a)), bringing in more than the minimum outdoor air-flow rate ( m oa ) required will increase the sensible and latentloads on the cooling coil. So MZHC decides to keep m oa atthe minimum; see Figure 10(d). Figure 12 shows the simulation results for a mild week.It is found that using
MZHC leads to ∼
60% energy savingswhen compared to BL , as shown in Figure 9. This significantreduction in energy consumption can be attributed to threemain reasons. Two of them are the same as those explainedin Section 6.1; the effects here are more prominent and thedetails are discussed in the subsequent paragraph. The third12
24 48 72 96 120 144 1681020300 24 48 72 96 120 144 1680501000 24 48 72 96 120 144 168
Time (Hours) (a) Outdoor weather data used in simulations (outdoor air temperature, out-door air relative humidity, and solar irradiance).
MZHCBL
24 48 72 96 120 144 168015300 24 48 72 96 120 144 168
Time (Hours) (b) Comparing the power consumptions (fan, cooling, and total reheatpower).(c) Aggregate conditions of floors/meta-zones 1, 2, and 3 (temperatures andrelative humidities) when using
MZHC and BL . The black dashed lines arethe thermal comfort limits. MZHCBL
Time (Hours) (d) Conditions at the AHU (supply airflow rate, outdoor airflow rate, con-ditioned air temperature, and conditioned air humidity ratio).
Fig. 12: Comparison of the two controllers for a mild week (Feb/19 to Feb/26, Gainesville, Florida, USA).is that, since the outdoor weather is mild and dry (see Fig-ure 12(a)),
MZHC also decides to use “free” cooling whenpossible by bringing in more than the minimum outdoor airrequired which leads to further reduction in cooling energyconsumption.Figure 12(c) shows the aggregate zone temperature( T z , f ) and relative humidity ( RH z , f ) for all the threefloors/meta-zones. Unlike the results for the hot-humidweek, the humidity constraint is found to be only intermit-tently active as the outdoor weather is relatively dry. Thisprovides more room for optimizing the conditioned air tem-perature, which has two important implications: (i) reduc-tion in the cooling energy consumption, and (ii) minimal re-heat energy consumption. For example, see 13:00-22:00 hin Figure 12(d) where T ca is at its higher limit, and also seeFigure 12(b) where P cc is significantly reduced and P reheat isalmost zero. The energy savings when using
MZHC is found to besignificant as can be seen in Figure 9. Since the outdoorweather is cold and dry, there is a lot of room for optimiz-ing the control commands. The reasons for energy reductionwhen using
MZHC are the same as those explained in Sec-tion 6.2, therefore we do not discuss them in detail here.
Comment . Recall that as mentioned in Section 4, the In-novation Hub building uses
Single Maximum algorithm forzone climate control as opposed to the
Dual Maximum algo-rithm presented here [22]. In the case of
Single Maximum al-gorithm, the minimum allowed value for the supply airflowrate has to be high enough so that the heating load can bemet with low enough supply air temperature to prevent ther-mal 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 strat-ification. Therefore using
Single Maximum algorithm leads13o higher fan, cooling, and reheating energy consumptions.This also implies that the energy savings by the proposedcontroller will be even higher.
The proposed control architecture is designed to addressa number of limitations in the existing literature on multi-zone building control using MPC. The main one is the re-liance on a high-resolution multi-zone model, which can bechallenging to obtain. A low-resolution model of the build-ing is more convenient since such a model can be identifiedin a tractable manner from measurements. The challengethen is to convert the MPC decisions, which are computedfor the fictitious zones in the model, to the commands for theVAV boxes of the actual building. The proposed architecturedoes that by posing this conversion as a projection problemthat uses not only what the MPC computes but also feedbackfrom zones. The result is a principled method of computingVAV box commands that are consistent with the optimal de-cisions made by the MPC without needing dynamic modelsof indiviudal zones. At the same time, the use of feedbackfrom the zones ensures that zone climate states are also closeto the aggregate climate states computed by the MPC.The positive simulation results provide confidence onthe effectiveness of the proposed controller especially be-cause of the large plant model mismatch. The simulationtestbed mimics a real building closely, including the het-erogenous nature of the zones in the building. Another ob-servation from the simulations that should be emphasized isthe need to incorporate humidity and latent heat. The indoorhumidity constraint is seen to be active when using the pro-posed controller, especially during hot-humid weather. With-out humidity being explicitly considered, the controller islikely to have caused high space humidity in an effort to re-duce energy use.There are many avenues for further exploration, such asexperimental verification, modifying the formulation to min-imize cost instead of energy and to include demand charges,improving methods to forecast disturbances, etc.
Acknowledgements
This research reported here has been partially supportedby the NSF through award
References [1] Serale, G., Fiorentini, M., Capozzoli, A., Bernardini,D., and Bemporad, A., 2018. “Model predictive control(MPC) for enhancing building and HVAC system en-ergy efficiency: Problem formulation, applications andopportunities”.
Energies, (3), p. 631. [2] Shaikh, P. H., Nor, N. B. M., Nallagownden, P., Elam-vazuthi, I., and Ibrahim, T., 2014. “A review on opti-mized control systems for building energy and comfortmanagement of smart sustainable buildings”. Renew-able and Sustainable Energy Reviews, , pp. 409 –429.[3] Goyal, S., and Barooah, P., 2013. “Energy-efficientcontrol of an air handling unit for a single-zone VAVsystem”. In IEEE Conference on Decision and Con-trol, pp. 4796 – 4801.[4] Joe, J., and Karava, P., 2019. “A model predictivecontrol strategy to optimize the performance of radiantfloor heating and cooling systems in office buildings”. Applied Energy, , pp. 65 – 77.[5] Raman, N., Devaprasad, K., Chen, B., Ingley, H. A.,and Barooah, P., 2020. “Model predictive control forenergy-efficient HVAC operation with humidity and la-tent heat considerations”.
Applied Energy, , Decem-ber, p. 115765.[6] Chen, X., Wang, Q., and Srebric, J., 2016. “Occupantfeedback based model predictive control for thermalcomfort and energy optimization: A chamber experi-mental evaluation”.
Applied Energy, , pp. 341 –351.[7] Ma, J., Qin, J., Salsbury, T., and Xu, P., 2012. “De-mand reduction in building energy systems based oneconomic model predictive control”.
Chemical Engi-neering Science, (1), pp. 92 – 100. Dynamics, Con-trol and Optimization of Energy Systems.[8] Bengea, S. C., Kelman, A. D., Borrelli, F., Taylor,R., and Narayanan, S., 2013. “Implementation ofmodel predictive control for an HVAC system in a mid-size commercial building”. HVAC&R Research, ,pp. 121–135.[9] Radhakrishnan, N., Su, Y., Su, R., and Poolla, K., 2016.“Token based scheduling for energy management inbuilding HVAC systems”. Applied Energy, , pp. 67– 79.[10] Png, E., Srinivasan, S., Bekiroglu, K., Chaoyang, J.,Su, R., and Poolla, K., 2019. “An internet of thingsupgrade for smart and scalable heating, ventilation andair-conditioning control in commercial buildings”.
Ap-plied Energy, , pp. 408 – 424.[11] Ma, Y., Richter, S., and Borrelli, F., 2012. “Chap-ter 14: Distributed model predictive control for build-ing temperature regulation”.
In Control and Opti-mization with Differential-Algebraic Constraints, ,March, pp. 293–314.[12] Mei, J., and Xia, X., 2018. “Multi-zone buildingtemperature control and energy efficiency using au-tonomous hierarchical control strategy”. In 2018 IEEE14th International Conference on Control and Automa-tion (ICCA), pp. 884–889.[13] Patel, N. R., Risbeck, M. J., Rawlings, J. B., Wenzel,M. J., and Turney, R. D., 2016. “Distributed economicmodel predictive control for large-scale building tem-perature regulation”. In 2016 American Control Con-ference (ACC), pp. 895–900.1414] Yang, Y., Hu, G., and Spanos, C. J., 2020. “HVACenergy cost optimization for a multizone building viaa decentralized approach”. IEEE Transactions on Au-tomation Science and Engineering, (4), pp. 1950–1960.[15] Yushen Long, Shuai Liu, Lihua Xie, and Johansson,K. H., 2016. “A hierarchical distributed MPC forHVAC systems”. In 2016 American Control Confer-ence (ACC), pp. 2385–2390.[16] Jorissen, F., and Helsen, L., 2016. “Towards an au-tomated tool chain for MPC in multi-zone buildings”.In 4th International Conference on High PerformanceBuildings, pp. 1–10.[17] Li, X., and Wen, J., 2014. “Review of building energymodeling for control and operation”. Renewable andSustainable Energy Review, , p. 517–537.[18] Kim, D., Cai, J., Ariyur, K. B., and Braun, J. E.,2016. “System identification for building thermal sys-tems under the presence of unmeasured disturbances inclosed loop operation: Lumped disturbance modelingapproach”. Building and Environment, , pp. 169 –180.[19] Zeng, T., and Barooah, P., 2020. “Identification of net-work dynamics and disturbance for a multi-zone build-ing”.
IEEE Transactions on Control Systems Technol-ogy, , August, pp. 2061 – 2068.[20] Guo, Z., Coffman, A. R., Munk, J., Im, P., Kuru-ganti, T., and Barooah, P., 2020. “Aggregation anddata driven identification of building thermal dynamicmodel and unmeasured disturbance”. Energy andBuildings , September. in press, available online Sept2020.[21] Fritzson, P., and Engelson, V., 1998. “Modelica —A unified object-oriented language for system model-ing and simulation”. In ECOOP’98 — Object-OrientedProgramming, E. Jul, ed., Springer Berlin Heidelberg,pp. 67–90.[22] ASHRAE, 2011. The ASHRAE handbook : Applica-tions (SI Edition).[23] Liang, W., Quinte, R., Jia, X., and Sun, J.-Q., 2015.“MPC control for improving energy efficiency of abuilding air handler for multi-zone vavs”.
Building andEnvironment, , pp. 256 – 268.[24] Taylor, S. T., 2007. “VAV system static pressure set-point reset”. ASHRAE journal, .[25] Jorissen, F., Reynders, G., Baetens, R., Picard, D., Sae-lens, D., and Helsen, L., 2018. “Implementation andVerification of the IDEAS Building Energy SimulationLibrary”. Journal of Building Performance Simulation, , pp. 669–688.[26] Roulet, C.-A., Heidt, F., Foradini, F., and Pibiri, M.-C., 2001. “Real heat recovery with air handling units”. Energy and Buildings, (5), pp. 495 – 502.[27] ASHRAE, 2017. The ASHRAE handbook fundamen-tals (SI Edition).[28] Goyal, S., and Barooah, P., 2012. “A method for model-reduction of non-linear thermal dynamics of multi-zonebuildings”. Energy and Buildings, , April, pp. 332– 340.[29] ASHRAE, 2016. ANSI/ASHRAE standard 62.1-2016,ventilation for acceptable air quality.[30] Goyal, S., Ingley, H., and Barooah, P., 2013.“Occupancy-based zone climate control for energy effi-cient buildings: Complexity vs. performance”. AppliedEnergy, , June, pp. 209–221.[31] National Solar Radiation Database (NSRDB). https://nsrdb.nrel.gov .[32] Andersson, J. A. E., Gillis, J., Horn, G., Rawlings,J. B., and Diehl, M., 2019. “Casadi: a software frame-work for nonlinear optimization and optimal control”.
Mathematical Programming Computation, (1), Mar,pp. 1–36.[33] W¨achter, A., and Biegler, L. T., 2006. “On the im-plementation of an interior-point filter line-search algo-rithm for large-scale nonlinear programming”. Mathe-matical Programming, (1), Mar, pp. 25–57.[34] Williams, J., 2013. Why is the supply air temperature55F? http://8760engineeringblog.blogspot.com/2013/02/why-is-supply-air-temperature-55f.htmlhttp://8760engineeringblog.blogspot.com/2013/02/why-is-supply-air-temperature-55f.html