Building Efficiency and Sustainability in the Tropics ( SinBerBEST ) Title Handling model uncertainty in model predictive control for energy efficient buildings Permalink

Model uncertainty is a significant challenge to more widespread use of model predictive controllers (MPC) for optimizing building energy consumption. This paper presents two methodologies to handle model uncertainty for building MPC. First, we propose a modeling framework for online estimation of states and unknown parameters leading to a parameter-adaptive building (PAB) model. Second, we propose a robust model predictive control (RMPC) formulation to make a building controller robust to model uncertainties. The results from these two approaches are compared with those from a nominal MPC and a common odeling uncertainty uilding control odel predictive control uilding HVAC system alman filter building rule based control (RBC). The results are then used to develop a methodology for selecting a controller type (i.e. RMPC, MPC, or RBC) as a function of building model uncertainty. RMPC is found to be the superior controller for the cases with an intermediate level of model uncertainty (30–67%), while the nominal MPC is preferred for the cases with a low level of model uncertainty (0–30%). Further, a common RBC outperforms MPC or RMPC if the model uncertainty goes beyond a certain threshold (e.g. 67%). © 2014 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license


Introduction
Reducing the energy consumption of buildings by designing smart controllers for operating the HVAC system in a more efficient way is critically important to address energy and environmental concerns [1]. Advanced control algorithms are considered major enablers to achieve higher energy efficiency in commercial buildings. Entire sections of the ASHRAE 90.1 standard [2] are dedicated to the specification of control requirements. Although the optimal control of an HVAC system is a complex multi-variable problem, it is standard practice to rely on simple control strategies that include on-off controllers with hysteresis, and PID controllers.
For optimal control design a thermal model of the building is needed. To achieve building-level energy-optimality, building model should be able to capture the interaction between physically connected spaces in the building, heat storage in walls, and provide an accurate prediction of temperature in the building. Control algorithm on the other hand, should be able to minimize energy consumption and optimize thermal comfort by exploiting occupancy schedules, weather forecast, and system dynamics (i.e. a model to predict temperature evolution of indoor air), and satisfy state (i.e. room air temperature and wall temperatures) and inputs (i.e. discharge air temperature and air mass flow rate) constraints and operate the HVAC system of the building in an optimal fashion within the range of operation of the components.
Model predictive control (MPC) is a promising control strategy that is capable of addressing all the aforementioned criteria and has shown results for achieving higher energy efficiency in buildings [3][4][5][6][7]. MPC can provide a potential building energy saving of 16-41% compared to the commonly used rule-based building HVAC controllers [3,8,9]. Other advantages of MPC for building HVAC systems include robustness, tunability, and flexibility [3]. Application of MPC for building energy control has been reported in the literature [3,8,7,[10][11][12][13][14]. There are different variations of nominal MPC such as distributed [13,15], robust [10,16] and stochastic [7,4]. MPC strategies to systematically address various challenges in building energy control. In [17] the authors propose a computationally tractable approximation of the nonlinear optimal control problem by which they optimize the predicted mean vote (PMV) index, as opposed to the static temperature range. A robust control strategy based on static pressure and supply air temperature reset control is presented in [18] for variable air volume (VAV) system. [19] proposes a controller based on a three mode robust control strategy where each mode addresses different control objectives and conditions; this proposed controller is robust in different load conditions. Authors of [20] showed that in presence of model uncertainty an H ∞ -robust controller achieves not only a robust performance on set-point tracking of the air-handling unit but also less energy consumption compared to the pole-placement controller. Authors in [21] observed that indoor zone volume acts as system's bifurcation parameter. A multi-variable regulation strategy based on feedback linearization is used to prevent secondary Hopf bifurcation. The designed control improves the limit cycle behavior and decreases indoor temperature variation.
However, these control techniques rely heavily on a perfect (or almost perfect) mathematical model of the building and a perfect estimation of the unmodeled dynamics of the system [3] to achieve considerable energy saving. In [8] the authors argue that, based on industrial experience, modeling is the most time-demanding and costly part of the automation process. Recently, numerous mathematical models of building thermal dynamics have been proposed in the literature. Resistor-capacitor (RC) models with disturbances to capture unmodeled dynamics have been proposed in [9,3,22]. A bilinear version of an RC model is presented in [7] that takes into account weather predictions to increase building energy efficiency. In [23], the authors found that time varying properties such as occupancy can significantly change the dynamic thermal model and influence how building models are identified. While modeling a multi-zone building, the authors of [23] observed that the experimental data often did not have sufficient quality for system identification and hence, proposed a closed-loop architecture for active system identification using prediction-error identification method (PEM). Although a great deal of progress has been made in modeling the thermal behavior of building envelope and HVAC system [9,22,3,24,5,6], the random nature of some components of these systems makes it very hard to predict, with high fidelity, the temperature evolution of the building using mathematical models. Buildings are dynamical systems with uncertain and timevarying physical and occupancy characteristics. The heat transfer characteristics of a building are highly dependent on the ambient conditions. For instance, heat transfer properties such as convective heat transfer coefficient h, of peripheral walls is dependent on outside temperature, wind speed and direction. Also, unmodeled dynamics of a building [3] is function of (1) external factors: ambient weather conditions such as radiative heat flux into the walls and windows, and cloudiness of the sky, and (2) internal factors: such as occupancy level, internal heat generation from lighting, and computers. These quantities are highly time-varying and therefore the dynamics of the building and, consequently, parameters of the mathematical model need to constantly adapt to this change over time.
One approach to increase the accuracy of the linear building models is to use an adaptive parameter estimation technique such that the building parameters are updated as the environment changes which leads to an adaptive modeling framework. Although this technique has been used for joint state-parameter estimation in other applications [25][26][27], to the best of the authors' knowledge, this paper is the first study on developing adaptive modeling framework for simultaneous estimation of building parameters, states and unmodeled dynamics.
Four approaches can be taken to model dynamic behavior of buildings and overcome model uncertainty for building controls: 1. Develop detailed nonlinear physical models for building [28,29] The first approach is typically computationally expensive. Consequently, its application for real-time building controls is limited.
The second approach provides accurate information about timevariation of influential factors on building performance but this approach is not cost-efficient and can be limited by the possibility of adding new sensors to a building. The third and fourth approaches are promising and they are the focus of this paper. In particular, we develop a parameter adaptive building (PAB) model and design a robust MPC for buildings. In this paper we build upon our previous work reported in [3,9,33,22,16].
The overall contribution of this paper is putting together modeling, control and co-design in a coherent framework to develop a methodology for selecting a controller type (i.e. RMPC, MPC, or RBC) as a function of building model uncertainty. Particular contributions are: The rest of the paper is organized as follows. Section 2 explains the experimental setup used to collect data for this study. We present the proposed parameter adaptive building (PAB) model and the developed parameter/state estimation technique in Section 3. Controller design and performance results for MPC and RMPC, as well as the indices based on which we assess the performances of the introduced controllers are presented in Section 4. Conclusions are drawn in Section 5.

Test-bed and historical data
The model studied in this paper is a model for an office room in Lakeshore building at Michigan Technological University. This room is surrounded by two rooms and a corridor in the building and connected to the outdoor area with a thick concrete wall and two south-oriented double-layered windows. Each room is equipped with temperature and humidity sensors (Uni-curve Type II) with the temperature accuracy of ± 0.2 • C as part of the building management system (BMS). We have used a different sensing device, (temperature data logger with accuracy of ±0.8 • C) to account for spatial temperature variation in the room and sensing accuracies of individual sensors. Location of the zone sensors are shown in Fig. 1. Temperature readings from these two sensors are shown in Fig. 2. We follow the methodology proposed in [33] to find the temperature measurement accuracy, which is obtained to be ±0.8 • C, and is used in the state estimation algorithm which is described in Section 3.3. Outdoor temperature is also measured by the BMS system.
The HVAC system in the building uses ground-source heatpumps (GSHP) to obtain required energy for heating purposes. Each Fig. 1. Location of the temperature sensors in the test-bed. The sensor 1-a is the room temperature sensor and the sensor 1-b is a temperature data logger installed to calculate measurement errors. unit in this system provides heating for an individual zone. Therefore, a unit operates when heating is required for its zone: the setpoint can be defined independently based on the functionality of each zone. The HVAC system uses an on-off controller to provide a desired temperature for each zone. Zone temperatures are measured with a sampling period of 60 s.

Parameter adaptive building (PAB) model
Building models proposed in the literature depend on many parameters. The reason is that buildings are composed of many sub-systems and a variety of thermal mechanisms takes place in the building such as heat conduction through walls, forced convection due to air conditioning systems, and thermal radiation from outside. A mathematical model that is descriptive enough to accurately explain these phenomena will end up with many time-varying parameters. Finding the best parameters at each time step is shown to be cumbersome [23]. In this section we propose and develop a novel parameter adaptive building (PAB) model that facilitates this parameter tuning process in an online and automatic fashion. The architecture of the proposed PAB model is shown in Fig. 3. Measurement data from various sensors such as temperature and airflow are stored in a data repository. The PAB model has a parameter update module which takes care of automatic parameter tuning on the fly, and is explained in detail later in this section. The PAB model works as follows: historical data is used to perform off-line, one-step model calibration. The obtained parameters from model calibration is used in the parameter update module (exploiting Kalman filtering algorithm) as an initial set of Here we first review fundamental heat transfer mechanisms in buildings, leading to a mathematical model of building climate, on top of which we develop the PAB model in the rest of this section. Fig. 4 depicts the schematic of a typical room studied in this paper. We use lumped model analysis [34] to reduce the complexity of the model, and obtain a low order model, suitable for control purposes. As a simplifying assumption, temperature is considered uniform inside the room. We use RC model from [22] in which the building is considered as a network of nodes. We account for time varying parameters by updating the parameters on the fly. More details on online parameter estimation is presented in Section 3.2.

Heat transfer
There are two types of nodes in the building network: walls and rooms. Consider in total n nodes, m of which represent rooms and the remaining n − m nodes represent walls. We denote the temperature of room i with T r i . The wall node and temperature of the wall between room i and j are denoted by (i, j) and T w i,j , respectively, thermal dynamics of which is governed by the following equation: where C w i,j , ˛i ,j and A w i,j are heat capacity, radiative heat absorption coefficient and area of wall between room i and j, respectively. R i,j k is the total thermal resistance between the centerline of wall (i, j) and the side of the wall where node k is located. Q rad i,j is the radiative heat flux density on wall (i, j). N w i,j is the set of all of neighboring nodes to node w i,j . r i,j is wall identifier which is equal to 0 for internal walls, and equal to 1 for peripheral walls (i.e. either i or j is an outside node). In Eq. (1) the left term denotes the rate of change of stored heat in the wall between room i and room j. The first term of the right hand side of this equation represents the flow of heat between room k and wall (i, j) due to temperature difference and the second term shows the heat flow to the wall, due to solar radiation. Temperature dynamics of the ith room is modeled by the following equation: where T r i , C r i andṁ r i are the temperature, heat capacity and air mass flow into room i, respectively. c a is the specific heat capacity of air, and T s i is the temperature of the supply air to room i. i is window identifier which is equal to 0 if none of the walls surrounding room i have a window, and is equal to 1 if at least one of them has a window. w i is the transmissivity of glass of window i, A win i is the total area of windows on walls surrounding room i, Q rad i is the radiative heat flux density per unit area radiated to room i, andQ int i is the internal heat generation in room i. N r i is the set of all of the neighboring room nodes to room i. In Eq. (2) the left term denotes the rate of change of stored heat in the air in room i. The first term of the right hand side of this equation represents the flow of heat between node k and room i due to temperature difference, the second term shows the heat flow delivered by the heating system, the third term represents the total radiative heat passing through the windows and the fourth term is the internal heat generation inside room i. More details of building thermal modeling and estimation of the unmodeled dynamics is available in [9,22,3]. Note that we approximate the values of Q rad i (t) andQ int (t) based on the following equations: where T out and « are the outside air temperature and CO 2 concentration in the room, respectively [35]. Air ventilation is considered constant as a simplifying assumption. A more sophisticated model for gas transport process in buildings can be found in [36]. Parameters , , and are obtained by the parameter estimation algorithm detailed in Section 3.3.
We model the radiative heat transfer between building and ambient environment as proposed in [37]. The amount of heat transferred from the building to the environment is given by the Stefan-Boltzmann law: where T bldg is the average temperature of the building. We also consider solar radiation heat transfer, Q solar absorbed by the walls, and the room through the windows. The data used in this paper is based on the past 30 years monthly average of solar radiation for flat-plate collectors facing south (resembling the south facing flat vertical walls of the building), and is obtained from NREL (National Renewable Energy Laboratory) [38] database for Houghton, MI in January. Furthermore, we take into account the radiation cooling at night (i.e. sky thermal radiation to the building) based on the proposed relation in [37]: where K is the coefficient related to the cloud height and C is a function of cloud coverage. We use K = 0.34 and C = 0.8 for simulations, based on the results in [37]. T out is the outside air temperature, and RH is the air relative humidity percentage. The total radiation exchange between building and ambient environment is then given by: Note that Q sky and Q solar are heat flow into the building, and Q bldg , is the heat flow from the building to the environment.

System dynamics
Heat transfer equations for walls and rooms yield the following system dynamics: where x t ∈ R n is the state vector representing the temperature of the nodes in the thermal network, u t ∈ R lm is the input vector representing the air mass flow rate and discharge air temperature of conditioned air into each thermal zone, and y t ∈ R m is the output vector of the system which represents the temperature of the thermal zones. l is the number of inputs to each thermal zone (e.g., two for air mass flow and supply air temperature). C is a matrix of proper dimension and the disturbance vector is given

Disturbance
Following the intuitive linear relation between outside temperature T out , internal heat generationQ int , and solar radiation Q rad , with the building internal temperature rise we approximate d t with an affine function of these quantities, leading to: where a, b, c, e are constants to be estimated. By substituting (3) and (4) into (9) and rearranging the terms, we get: where a = a + c, b = b , and e = a + b + e. Therefore, only measurements of outside air temperature and CO 2 concentration levels are needed to determine the disturbance to the model. The values of a, b, and e are estimated along with other parameters of the model.

Additive uncertainty
We linearize the original nonlinear dynamic system and use Euler's discretization method to obtain a linear discrete-time system. We also add an additive uncertainty to the state update equation to account for model uncertainties, leading to: where the uncertainty w k ∈ R r is a stochastic additive disturbance. t ∈ R refers to time in continuous-time domain and k ∈ Z refers to time in discrete-time domain. The set of possible disturbance uncertainties is denoted by W and w k ∈ W ∀k = 0, 1, . . ., N − 1.
For this study, we consider box-constrained disturbance uncertainties given by

State-parameter estimation
Using (1) for each wall and (2) for each room node in the building network, system dynamics is given by: where x 1 is the room temperature (T r 1 ), and x 2 , x 3 , x 4 , x 5 are the peripheral walls' temperature (i.e. T w12 , T w13 , T w14 , T w15 ). T 2 , T 3 , T 4 , T 5 are the temperatures of the surrounding zones, as shown in Figs. 4 and 5. These temperatures act as disturbance to the system dynamics for a single zone thermal model, and x is the state vector: One way to adapt the model to account for time varying parameters is to assume that all the parameters of the model are independent, and hence define a state corresponding to each parameter. However, this would lead to excessive number of states (e.g. 18 states for a room shown in Fig. 4). To overcome this problem, we take a different approach. We reduce the number of states by exploiting the redundancies in the resulting model. For instance, thermal properties of wall material (e.g. specific heat capacity and conductive heat transfer coefficient) are the same across the building, as these are functions of the materials used as the building walls. In addition, the thickness of internal walls and thickness of peripheral walls are the same throughout the building. Following this approach, we are able to reduce the number of independent parameters from 18 to 10. Hence we re-write the thermal equations of the walls, i.e. (13b)-(13e) as follows: As shown in (19), C w R w is not a function of the area of wall: where c w , k w , A w and L w are the specific heat capacity, conductive heat transfer coefficient of wall material, area and thickness of wall, respectively, and h in is the indoor convective heat transfer coefficient. Hence, we can use one common term to express thermal capacitance-resistance between centerline of each wall and the node on each side of the wall for the equations of walls in the building.
We designate a state variable to all the independent timevarying parameters of the system as follows: Rate of change of these states is equal to zero, as shown in the corresponding state update Eq. (30). We then add a low-magnitude fictitious noise to the dynamics of parameters to allow slow changes in their values over time.
u is the input vector given by: In summary, we express the dynamics of the system using following state update model: where w k and v k are the process and measurement noise and are assumed to be zero mean multivariate Gaussian process with variance W k and V k , (i.e. w k ∼N(0, W k ) and v k ∼N(0, V k )), respectively.

Estimation algorithm
In order to estimate the unknown parameters of the system we augment the states of the system with a vector p k which stores the parameters of the system, with a time evolution dynamics of p k+1 = p k , as will be detailed in Appendix B. Due to the multiplication of states and parameters the resulting dynamic system is nonlinear. Nonlinear estimation algorithms such as extended Kalman filtering (EKF) or unscented Kalman filtering (UKF) can then be exploited to simultaneously estimate the states and the parameters of the system.
An alternative to using a Kalman filter would be a simple observer. However, given the random variations, inaccuracies and uncertainties in the system dynamics, as described earlier in the paper, using a Kalman filter is suggested in order to get a statistically optimal estimate of system states [39,40].
In our previous work [41] we showed that UKF outperforms EKF for building parameter estimation. Thus, we only focus on UKF in this study. We present an algorithmic description of the UKF in Appendix B, omitting some theoretical considerations.

Estimation results
The test-bed from Section 2 was used to collect measurements from January 11 to January 24, 2013. To remove noise from the temperature measurements, a second order Butterworth lowpass filter with cutoff frequency of 0.001 Hz was used. Fig. 5 shows the temperatures of the neighboring zones and the outside temperature which act as disturbance to the PAB model. Fig. 6 depicts the model inputs including the air mass flow rate and the supply air temperature. In order to obtain the best initial parameter values for the Kalman filter algorithm, we first perform a (static) parameter identification on the historical data. We consider the first part of the data as training set (shown in red in Fig. 7), and obtain the best parameters that minimize the least square error between the simulation and the measurement data. The result of this step is used to simulate the temperature evolution of the room air for the next three days (shown in black in Fig. 7). Due to time-varying parameters and disturbance to the model, it is difficult to find a set of parameters for the model which results in good temperature tracking for all days including weekdays and weekends, and hence, as shown in Fig. 7, the results of simulations for the following days in the testing data set is even worse.
The obtained initial parameters from the off-line calibration step is used as initial value for the UKF algorithm. For the off-line parameter calibration practice, we used the historical data of two weeks where the first 60% of the data was used for training (calibration) and the remaining 40% of data was used for testing. The due to the changing environment. Note that the first part of the estimation of wall temperature by UKF leads to overshoot in the wall temperature, however, this overshoot is quickly recovered as UKF uses more data to tune the parameters more accurately.
UKF is also tested to estimate the temperature in the presence of process and measurement noise (w and v, respectively) as shown in Fig. 11. We add process and measurement noise to the model and use UKF to estimate the temperatures. UKF is used to estimate the temperature from the measurements. Performance of UKF is shown with model uncertainty w, and measurement noise v, given by w k ∼N(0, 0.2) and v k ∼N(0, 1.4), respectively. As seen in Fig. 11, UKF is able to cancel out the effect of noise very effectively.

Controller design
In this section we study the impact of the use of the PAB model in a model-based control design framework. State-of-the-art is to use a fixed-parameter model to design MPC for buildings. We propose using the updated parameter model obtained using the Kalman filter estimation process at each time step as shown in Fig. 3, which results in a more accurate model and hence lower model uncertainty. The underlying assumption here is that the parameters of the system do not change from time t to t + 1. At the next time step, MPC uses the model with updated parameters, to derive the optimal inputs. Inputs are implemented on the system and at the next sampling time new states (temperatures) are measured and sent to the PAB model, and this process repeats. We also formulate a nominal MPC and a robust model predictive control (RMPC), and study their performances for various model uncertainty levels. MPC assumes that the model is perfect (no uncertainty), and the RMPC assumes that the model is uncertain and designs a robust control policy for a specific class of uncertainty. The results from MPC and RMPC are compared to a conventional rule-based control (RBC) for a typical building. Novel performance indices are proposed to compare the performance of these controllers. We also present a methodology to select the best controller among the ones studied in this section for any given model uncertainty, which leads to optimum trade-off between energy consumption and comfort level.

ASHRAE requirements for building climate control
ASHRAE's Standard 55 [42], Thermal Environmental Conditions for Human Occupancy, suggests the condition which is acceptable to at least 80% of occupants. According to this standard, the ideal temperature in typical clothing in summer (0.35-0.6 clo) is in the range of 22.5-26 • C. The operative temperature for occupants in normal clothing insulation in winter which is between 0.8 and 1.2 clo should be in range of 20-23.5 • C. This temperature range is based on a metabolic rate of 1.2 met (70 W/m 2 ) and 60% RH. More details can be found in [43][44][45]. ASHRAE's Standard 62.1 [46], Ventilation for Acceptable Indoor Air Quality, explains outdoor air ventilation requirements for different types of indoor spaces. When the major contamination source is proportional to number of occupants, the minimum ventilation rate is enforced in CFM (L/s) and when other factors play the main role in contamination, the minimum ventilation rate is enforced in CFM/ft 2 (L/s m 2 ) [45]. We use this as a guideline for control design in this section.

Rule-based control (RBC)
The rule-based controller in this study is a conventional on-off HVAC controller. The time constant of the control action implementation is t = 1 h. The controller opens the dampers of conditioned air flow to the thermal zones when heating is required and keeps it fully open for the duration of t. In the next time step the controller checks the temperature again and adjusts the damper position if the room temperature is within the comfort zone, or keeps it open if the room air temperature is still outside the comfort zone. In on-off control, position of the dampers can be either the min value or the max value. When system goes to the cooling mode, supply air temperature changes accordingly. The experimental data presented here is for the heating mode only. To be consistent and to perform a fair comparison, we use the same time constants t for all controllers.

Model predictive control (MPC)
A model predictive control problem is formulated with the objective of minimizing a linear combination of the total and the peak airflow. We implement the control inputs obtained from the MPC with the linearized system dynamics of the model on the original nonlinear model for forward simulation.
Fan energy consumption is proportional to the cubic of the airflow. Hence minimizing the peak airflow would dramatically reduce fan energy consumption. We have considered a cost function for the MPC which comprises linear combination of the total input airflow ( 1 norm of input) and the peak of airflow ( ∞ norm of input). The alternative would be to use the actual nonlinear function of fan energy consumption. However, it would lead to nonlinear MPC which is much slower than linear MPC. We use the proposed cost function to achieve better computational properties. Also in order to guarantee feasibility (constraint satisfaction) at all times, we implement soft constraints. The predictive controller solves at each time step the following optimization problem:  y t+k|t = Cx t+k|t , k = 1, . . ., N (33d) where U t = [u t|t , u t+1|t , . . . m, u t+N−1|t ] is vector of control inputs, and = [ε t+1|t , . . .m, ε t+N|t ] and = [ε t+1|t , . . .m, ε t+N|t ] are the slack variables used to utilize soft constraints on room temperature. y t+k|t is the room temperature vector, d t+k|t is the disturbance load prediction, and T t+k|t and T t+k|t for k = 1, . . . m, N are the lower and upper limits on the room temperature, respectively. U t+k|t and U t+k|t are the lower and upper limits on the airflow input by the variable air volume (VAV) damper, respectively. Note that based on the ASHRAE Standard 62.1 -Section 6.2.6.1, during unoccupied hours, ventilation systems should be able to maintain the required non-zero ventilation rates (U t+k|t > 0) in the breathing zone [46]. is the penalty on the comfort constraint violations, and Ä is the penalty on peak power consumption.
At each time step only the first entry of U t is implemented on the model. At the next time step the prediction horizon N is shifted leading to a new optimization problem. The prediction horizon is N = 24, and at each time step only the first entry of the input vector U t is implemented on the model. This process is repeated over and over until the total time span of interest is covered. We use YALMIP [47] to set up the MPC problem in MATLAB.

Robust model predictive control (RMPC)
We consider additive uncertainty to the system model as previously described in (11). A schematic of the robust optimal control implementation on the nonlinear building model is shown in Fig. 12. In RMPC algorithm, the cost function is the same as in the one in MPC case: However, state and input constraints are as follows: x t+k+1|t = Ax t+k|t + Bu t+k|t + E(d t+k|t + w t+k|t ) k = 0, 1, . . ., N − 1 (35a) y t+k|t = Cx t+k|t k = 1, 2, . . ., N (35b) The only difference with respect to MPC algorithm is the introduction of additive uncertainty term w in the state update equation.
Using this formulation, we derive a robust counterpart of an uncertain optimization problem in which constraints are satisfied for all possible uncertainties, and worst-case objective is calculated.
It is shown in [16] that the open-loop constrained robust optimal control (OL-CROC) is conservative. The closed-loop constrained robust optimal control (CL-CROC) formulation overcomes this issue but it can quickly lead to an intractable problem [48]. Next, we review the feedback prediction concept followed by our proposed formulation to improve upon the feedback prediction scheme.

Feedback predictions
The idea in feedback prediction, is to introduce new decision variables and parameterize the future control sequences using the future disturbances and an additive independent decision variable.
Define an affine disturbance feedback as:  (37) and the vector of disturbances is given by w = [w 0 w 1 · · · w N−1 ] .
The control sequence is parameterized directly in the uncertainty. What we have here is basically a sub-optimal version of the closed-loop min-max solution [48].

Two-lower-diagonal structure (TLDS)
One of the parameterizations introduced in [48] is Lower Triangular Structure (LTS). The main problem with the min-max formulation based on LTS parameterization is the excessive number of decision variables and constraints. The reason is the highdimensional parameterization of matrix M. To resolve the issue of high-dimensional parameterization of matrix M, we propose the following new parameterizations.
By analyzing the structure of the optimal matrix M, it was observed that the parameterization of the input does not need to consider feedback of more than past two values of w at each time, hence we propose the following disturbance feedback. and the corresponding parameterization matrix M is an N × N matrix that has the entries on the first and second diagonal of M below its main diagonal as decision variables and 0 elsewhere. n remains as in (37). With this structure we exploit the sparsity of the feedback gain matrix to enhance the computational characteristics of the controller.

Performance indices
To compare the overall performance of the proposed controllers we define indices to measure the energy consumption and comfort level provided by each controller. In addition, we define a new index to evaluate the overall performance of each controller considering both the energy and comfort indices.
• The energy index I e in (kWh) is defined: where cooling power P c , heating power P h and fan power P f are determined by: c p = 1.012 (kJ/kg • C) is the specific heat capacity of air and 3 = −6.06 × 10 −13 (CFM −3 ) and ˛2 = 0.73 × 10 −8 (CFM −2 ) and 1 = −1.2 × 10 −3 (CFM −1 ) and ˛0 = 59.2 and ˙v =ṁ/ is in (CFM), where is density of air and P f (t) is in % of nominal power consumption [49]. Using these constants, the fan power values, in (kW), can be calculated.
• The discomfort index I d in degree Celsius hour ( • Ch) is defined as the sum of all the temperature violations in the course of a day.
is the comfort zone at time t and 1 is the indicator function. • A good control performance means not only low energy consumption, but also low resulting discomfort. To assess the overall performance of the controllers, we need to examine both I e and I d at the same time. Using the two indices defined above we define a third index called overall performance index (I OP ). The intuition behind this new index is to take into account the energy and discomfort index in one single term. I OP is defined as: where I * d is the maximum allowed discomfort and || . || ∞ denotes infinity norm or the maximum value of energy indices among all three controllers. Negative value of I OP means that the discomfort index is not within the preferred range. The lower the I d and I e are, the higher the I OP will be. Therefore, the higher the I OP , the better the overall performance. In this study, the limit on the allowed discomfort index is heuristically chosen to be I * d = 0.5 ( • Ch) to ensure adequate comfort level.

Control results
To illustrate the effectiveness of the controllers proposed in where is the ∞ norm bound of the uncertainty as given by (12) and d = [d 1 , d 2 , . . ., d N ] is the disturbance realization vector. d represents transpose of vector d.
A time constant of t = 1 (h) is used for all controllers. We implement the introduced model predictive controllers with a prediction horizon of N = 24. The choice of N = 24 is to provide a good balance between performance and computational cost for the MPC framework in this study.
We use the following numerical values for parameters in Optimal controller and the resulting room temperature with the presence of a box-constrained uncertainty in four cases are depicted in Fig. 13. Measurements, as shown in black, shows the air mass flow and temperature recording for the room using a simple existing control policy of the building HVAC system. RBC represents the result of the rule-based control. MPC refers to the performance of a model-based control algorithm in which no knowledge of the model uncertainty is known a-prior to the control algorithm. RMPC refers to the simulation of the control algorithm which considers the model uncertainty bound and utilizes the uncertainty feedback strategy of (36) in designing the control policy.
We consider stochastic uncertainties with different uncertainty bounds ( ) as introduced in (12). The MPC does not have any a-priori information regarding the additive uncertainty, and calculates the controller solely based on the deterministic system dynamics. However the RMPC integrates the uncertainty bound information in the control derivation. Controller performances are evaluated based on indices introduced in Section 4.5. Problem is solved using CPLEX 12.2 [50] on a 2.67 GHz machine with 4 GB RAM.
Here are the discussions of the results:

Computational aspects
Exploiting the TLDS structure results in the same control law that was obtained from the LTS structure. However, matrix M of LTS has l . m . r((N(N − 1))/2) variables (quadratic in N) while matrix M of TLDS has l . m . r(2N − 3) variables (linear in N), and hence exhibits shorter computation time. On average, the simulation time for TLDS is 30% less than the LTS structure, as shown in Table 1.

Comfort
It is observed from Fig. 13, that the RMPC is the only controller that is able to keep the temperature within the allowed comfort zone, at all times during this test simulation, meaning maintaining minimum level of discomfort (I d ≤ I * d ), while RBC still does a very good job and MPC fails to do so, resulting to I d > I * d for all ı≥ 40 %. Fig. 14 depicts how discomfort index I d , varies with additive model uncertainty ı for MPC, and RMPC. Note that different data points for one ı value refers to simulations with different random sequences. The reason for such a wide variation of the simulation results, specially for large values of ı stems from the fact that depending on the value of the random variable at any time, the resulting disturbance vector can either lead to temperature rise or fall with respect to the nominal disturbance value. It is shown that RMPC manages to keep the perfect comfort level (I d = 0), for additive model uncertainty up to ı = 75%, while the MPC maintains the perfect comfort level for uncertainty bounds up to ı = 20%. The discomfort index for MPC goes as high as 4.61 • Ch while the value for RMPC reaches 1.2 • Ch in the worst case in the simulations corresponding to ı = 100%. Since RBC is not a model-based control technique, its performance does not depend on values of ı, hence the straight horizontal line in Fig. 14 Fig. 15 depicts the variations of energy index I e , versus the uncertainty bound on the unmodeled dynamics. It is clear that the energy index for RMPC increases dramatically with ı, while the energy index for MPC only changes slightly. However, this comes with the   Fig. 15 also shows energy consumption of RBC (I e = 1.43 × 10 4 kWh). MPC for all values of ı leads to a lower amount of energy consumption than RBC, but RMPC leads to more energy consumption than RBC soon after ı = 35%.

Comfort-energy trade-off
An important point to notice from Fig. 15 is how much more energy needs to be supplied to the HVAC system to maintain the comfort level in the presence of imperfect and faulty unmodeled dynamics predictions. Consider the case where ı = 75%. MPC will lead to a discomfort index of 1.7 • Ch on average, while the RMPC is able to maintain the temperature below a discomfort index of 0.016 • Ch on average. However this level of comfort provided by the RMPC comes at a cost of energy consumption of 3 times more than that of the MPC case. Note that due to the trade-off between comfort and energy consumption, the choice of which controller to use is on the building HVAC operator, and depends on various factors such as criticality of meeting the temperature constraints for the considered thermal zone in the building, and availability and price of energy at that time of the day/year. As observed from Figs. 14 and 15 the behavior of controllers vary considerably as the model uncertainty increases. For instance, the energy required to keep the same level of comfort for RMPC in the case of ı = 75% is almost 3 times the energy required to provide the same level of comfort when ı = 25%. Figs. 14 and 15 show the importance of a good model like PAB in minimizing the energy consumption of building HVAC systems for a desirable comfort level using model-based control techniques by accurately capturing the dynamics of the system. Fig. 16 demonstrates savings of MPC and RMPC versus RBC. As shown, the maximum theoretical energy saving of MPC compared to RBC is 36%, and that of RMPC is 30% for the building studied. These values decrease as model uncertainty increases. Energy saving of MPC versus RBC stays positive even for large values of model uncertainty, while energy saving of RMPC versus RBC is positive only for model uncertainty values up to about 34%, and is negative for larger model uncertainties (i.e. RMPC consumes more energy than RBC).

MPC and RMPC versus RBC
The results of an extensive study in [51] show that MPC HVAC control can potentially provide 16-41% building energy saving compared to rule-based controllers, which is in agreement with our findings. The saving depends on various factors including climate zone, insulation level, and construction type. Stochastic MPC was shown in [51] to be superior to the rule-based control given the uncertainties in occupancy and weather forecast. Our findings also show that the robust MPC outperforms the rule-based control in terms of energy consumption and user comfort. Although these two MPC techniques (robust and stochastic MPC) both address model uncertainty, they are formulated differently and hence can lead to different performance results. Given the accuracy of the PAB for removing model uncertainty, designing MPC scheme based on PAB is a promising solution for building control problem.
For simulation evaluation of energy consumption and provided comfort level, we have compared the overall performance of the three controllers using I OP . The results, as shown in Fig. 17, suggest that for model uncertainties less than 30% MPC is the best controller type. For model uncertainties above 30%, RMPC and RBC are close in performance while for ı between 30% and 67% RMPC is the best, and for model uncertainties larger than 67%, RBC leads to better overall performance than model-based control techniques. This information can be of utility for choosing a controller type for building HVAC system. As described in the paper, proper choice of building HVAC control would depend on the accuracy of the given building model. Range of uncertainties for a given building model can be obtained by taking the difference of the temperature predictions from the building model and temperature measurements from a building. The statistics of such uncertainty can be found once such data is available. The mean and variance of the uncertainty from the statistical analysis can be used to select the best controller type.

Summary and conclusion
Model uncertainty is an unavoidable challenge for modeling and model-based control of a building HVAC system. In this paper, we characterized the impact of model uncertainty on MPC controllers and presented two approaches to minimize model uncertainty for building controls. First, we presented a new modeling framework for simultaneous state estimation and parameter identification of building predictive models. This resulted in a parameter-adaptive building (PAB) model which captures system dynamics through an online estimation of time-varying parameters of a building model. The PAB model aims at reducing model uncertainty and can be used for both modeling and control. Second, we presented an MPC framework that is robust against additive uncertainty. The new framework is a closed-loop Robust Model Predictive Control (RMPC) utilizing uncertainty knowledge to enhance the nominal MPC. The RMPC is capable of maintaining the temperature within the comfort zone for model uncertainties up to 75%. The specific findings are listed below: 1. We constructed a nonlinear state space model by augmenting the parameters of the system into the state vector. We exploited the similarities in the physical properties such as wall materials and thicknesses in the building under study, and reduced the number of independent parameters in the building model. A similar approach is expected to apply to other building modeling practices. 2. We presented a PAB modeling framework that uses an unscented Kalman filter (UKF) to simultaneously estimate all the states of   heating air mass flow rate (kg/s) m r i air mass flow rate into the ith room (kg/s) N w i,j set of all of neighboring nodes to node w i,j N r i set of all of neighboring nodes to node i P c cooling power (kW) P f fan power (kW) P h heating power (kW) Q rad i,j radiative heat flux density on wall (i, j) (W/m 2 ) Q bldg heat flow from the building to the environment (W/m 2 ) Q sky heat flow from the sky to the building (W/m 2 ) Q solar solar heat flow to the building (W/m 2 ) Q int i internal heat generation in room i (W) r i,j wall identifier R i,j k total thermal resistance between centerline of wall i & j & the side of wall where node k is located (K/W) R win i,j thermal resistance of window between room i & j (K/W) RH relative humidity CO 2 concentration in the room (ppm) Stefan-Boltzmann law constant (W/m 2 K 4 ) w i transmissivity of glass of window i T bldg average temperature of the building ( • C) T r i temperature of the ith room ( • C) T out outside air temperature ( • C) T w i,j temperature of the wall between room i & j ( • C) T s i supply air temperature of the ith room ( • C) T lower limit on the room air temperature ( • C) T upper limit on the room air temperature ( • C) v measurement noise w i window identifier w k model uncertainty

Appendix B. Unscented Kalman filter
To perform UKF, we conduct the following initialization: x 0 = E (B.1) Each step of the UKF can be summarized as follows: Unscented Kalman filter algorithm Prediction: Calculate sigma points: Propagate each column of X k−1 through time: A-priori state estimate: A-posteriori state estimate: A-posteriori estimate of error covariance: P k = P − k − K k Pẑ kẑk K T k where: wherex − denotes a-priori estimate of state x. = (L + ), and = ˛2(L + ı) − L are the composite scaling parameters. ˛ is a scaling parameter that determines the spread of the sigma points aroundx, and is usually set to a small positive value (e.g. 1e − 4 ≤ ˛ ≤ 1). ı is a secondary scaling parameter which is usually set to 0 or 3 − L [40]. Q k is the process error covariance matrix and R k is the measurement noise covariance matrix. W where ˇ is a parameter used to incorporate the prior knowledge of the distribution of x. We use ˇ = 2 which is optimal for Gaussian distributions [52].