A Hybrid Physics-Based and Data Driven Approach to Optimal Control of Building Cooling/Heating Systems

This work integrates a physics-based model with a data driven time-series model to forecast and optimally manage building energy. Physical characterization of the building is partially captured by a collection of zonal energy balance equations with parameters estimated using a least squares estimation (LSE) technique and data initially generated from the EnergyPlus building model. A generalized Cochran-Orcutt estimation technique is adopted to describe the data generated from these simulations. The combined forecast model is then used in a model predictive control (MPC) framework to manage heating and cooling set points. This work is motivated by the practical limitations of simulation-based optimizations. Once the forecast model is established capturing sufficient statistical variability and physical behavior of the building, there will be no more need to run EnergyPlus in the optimization routine. The proposed methodology lends itself for real-life implementation of building energy management systems where predictive control is desired to reduce energy use and avoid demand charges and occupant discomfort. At each time step, it determines the optimal set point values of all building's zones and updates these values over time. In practice, the proposed control strategy can be implemented in commercial smart energy boxes to optimally control total daily energy-use costs.

statistical time series approach. In practice, the initial training of the model parameters is carried out using building simulation data. These estimates can then be refined and updated using actual data from the building of interest. The formulation of the MPC algorithm includes a multi-objective mathematical programming model that minimizes total operating energy cost and daily as-used demand charges as well as total deviation from thermal comfort bounds. Optimal control strategies (i.e., for the zonal set points) are updated every time that new estimates are available; this ensures that internal zone temperatures remain within prespecified limits over time. The novelty of this paper is on the specific combination and application of existing methods to optimize energy control of large buildings which are subject to stochastic externalities. In particular, the methodology integrates a physics-based zonal model with an advanced time series model to ensure enhanced accuracy and sensitivity of energy forecasts to incremental changes in control variables. Internal temperature measurements of different zones in the building are used in calculations. Initial training of the model parameters is carried out using highly granular building energy simulation (EnergyPlus or similar models). Unlike the current practices that run different scenarios to deal with stochastic externalities, the proposed forecasting model is adaptive and uses actual measurements to refine and update its forecast values. The basic idea of MPC is to form a model that is able to represent the future behavior of building cooling/heating dynamics and to provide optimal control actions for a specific time horizon [1]- [3]. Different modeling approaches can be employed for implementation of an MPC strategy. The first approach is based upon detailed physical modeling of cooling and heating dynamics. In this approach, physical characteristics of a building and its HVAC components are extracted and fed into a series of energy balance equations. The balance equations are then used for prediction of the future evolution of cooling/heating dynamics. For instance, references [4]- [7] present models for the cooling/heating dynamics using the resistance-capacitance (RC) network analogy. In these works, thermal capacitance and resistance between different components of a building or a zone, such as air, inside and outside surfaces of walls, windows and ceilings as well as heat flux due to solar radiation are represented through RC-network diagrams and heat transfer equations. The equations are then employed for prediction of cooling/heating dynamics and optimization of setpoint values. Models presented in [6]- [8] account for the effect of other dynamic variables or physical components such as occupancy, relative humidity, chilled/hot water temperature supplied to or returned from each building/zone, thermal storage equipment, etc.
The physical MPC approach is typically too complex to solve analytically for large granular building models. Therefore, researchers often apply physical MPC models for simpler problems, e.g., a single zone or a single room [3], [4]. An alternative is to use a data-driven approach to MPC by simply fitting a metamodel to the cooling/heating data regardless of the particular physical structure of the building. A number of studies have focused on linear statistical models where input and output data have linear forms. Autoregressive with exogenous variable (ARX) and autoregressive moving average with exogenous variable (ARMAX) are two examples of this approach [9], [10]. State-space modeling is another example of such an approach where model parameters are estimated over a number of specified system states [11], [12]. In other studies, researchers employ soft computing techniques, particularly artificial neural networks (ANNs), to address the complexity of building energy forecast modeling and optimization [13], [14]. ANNs provide linear or nonlinear model-free structures for prediction of energy demand using different input vectors, i.e., weather conditions or wall insulation thickness [15]. These predictions can then be used for optimization of cooling/heating loads [13].
Although data-driven models are typically simple to use, their implementation is often accompanied by a number of problems that may negatively influence their performance. ARX and ARMAX models, for instance, follow linear autoregressive structures that are not necessarily able to explain full variations of load dynamics. In addition, most soft computing techniques cannot guarantee full capture of complex interactions amongst building components, dynamic variables, and cooling/heating load, especially when the available real data is limited and the building includes a complicated multi-zone structure. To overcome such problems, a number of researchers employ a simulation-based approach to capture the dynamic behavior of buildings and thereby optimize energy use [16]- [18]. In this approach, first, a highly granular physics-based simulation model of the building is developed. Then, by designing and running different experiments, the behavior of building energy systems is captured over time. Simulation-based optimization techniques are applied to minimize the energy used to meet cooling/heating loads. However, this is costly and time consuming particularly for larger models where simulation optimization must run over a large spectrum of possible scenarios and run in near real time. Consequently, a number of recent works combine the benefits of both simulation and data-driven approaches to provide fast and effective control solutions [9]. These approaches are similar to that developed and used herein.

II. MODEL FRAMEWORK
Our proposed framework for optimal control of a building heating/cooling system is presented in Fig. 1. Model execution consists of two main phases: an offline phase and an online phase. The offline phase includes analysis using a set of historical data either generated by EnergyPlus, or directly gathered from buildings with the objective of constructing and training a heating/cooling dynamic model for the building. This heating/ cooling model is then used in a building energy forecast model to calculate the 24-h look-ahead forecast values for building energy consumption. In the online phase, both the heating/cooling model and energy forecast model are fed into an MPC scheme that is designed to provide the optimal cooling/heating set points for the next 24 h ahead. The MPC scheme is based on a dynamic programming approach, which runs every time that it receives new actual data from the building under study.
The following algorithm summarizes the proposed steps: This is a set of heat balance equations, which provides explicit relationships between zonal internal temperatures and effective power rate. The effective power rate represents the amount of cooling/heating rate in kW that the HVAC system supplies to each building zone during any specific time period. The model is used to forecast the k-step-ahead internal temperature for each zone.

3) Create an energy forecast model:
This combines the model from Step 2 with a time series model, and returns 24-h look-ahead forecast values for building total energy demand. We apply a generalized form of Cochran-Orcutt technique to estimate the model parameters [22]. Once the total energy demands for the next 24 h are forecasted, then one can calculate the total energy operating costs, which are then used in the next step.

4) Develop and run an MPC-based optimal control strategy:
We formulate a multi-objective dynamic programming code to search for optimal control set points for the next 24-h. The total operating costs of Step 3 as well as the total penalty for exceeding the thermal comfort bounds define the objective functions, and the heat balances in Step 2 form the state constraints. In conjunction with the thermal comfort objective function, an additional set of constraints are imposed to maintain the thermal comfort between specified bounds. As shown in Fig. 1, we find, at time , the optimal control set point values for times . Then, in the next hour, we update the optimal solutions for , when we receive feedback of new information on the building (i.e., from the building energy management system or real building).

III. HEATING/COOLING MODEL
Assume that a day is divided into a set of discrete time periods, . Then, according to the first law of thermodynamics, the total energy exchange associated with the th thermal zone, , at time step is given by (1) where , and , are the amounts of input and output energy at step , and is the amount of energy gained or lost at time by Zone . In this study, we assume that , is a function of internal and external temperatures [23] and can be calculated as follows: (2) and are the internal and the external temperature of the th zone at time , respectively. Assuming that the th zone is air conditioned by an HVAC system with effective power rate of (in kW) at time , we can rewrite (1) as follows: where is a white noise representing the additional unpredictable effects of convective internal loads, convective heat transfer from the zone surfaces, interzone air mixing effects and occupancy. In this study, we assume that such additional effects are normally distributed. In addition, according to [23], the internal temperature of the th zone at time can be written as (4) where is the duration of time period and is set equal to 1 h in this study. Furthermore, and are the heat capacity and the mass of air of the th zone in J/kg C, and kg, respectively. The unit of is and is C. Combining (3) and (4), we have (5) are independently and identically distributed. and cannot easily and accurately be determined in practice, since the mass of air and heat capacity are different for multiple points in a zone. Rather, the thermal balance model presented in (5) can be explained in terms of lost and delivered energy and the internal and external temperatures of the th zone as follows: (6) where represents the amount of unit increase (decrease) in the th zone internal temperature by one unit increase (decrease) in effective power rate over a time period. Hence, the effects of and are hidden in and , which are explicit and can be estimated using statistical techniques. For cooling seasons, it is logical to assume that and for heating seasons . Equation (6) can then be used to forecast the th zone internal temperature at time given that and are known [23]. The forecast values are calculated by (7) where is the forecast of internal temperature for the th zone at time . Assume that is the set point value of the zone and that ; where and are the lower and upper values for the possible set poins. In optimization phase, we set and find the corresponding by (7). and can be estimated using the least squares estimation (LSE) technique by minimizing , where and are the vectors of actual and forecasted internal temperature values for the th zone, and is the -norm of x. and can be calculated using numerical methods (see, e.g., [23]) or solving the following equations: Let and and (8) and (9) to zero; then we have (10) letting and , then, (10) can be rewritten in matrix form (11) where is the 2 1 vector of coefficients and . Equation (11) is given by multiplying both sides of (10) by the inverse of . According to (11), and and , are obtained as follows: (12) As mentioned in the previous section, the effective power rate, supplied to the th zone is not often measurable directly in real-world applications. We note that can be different from the electrical power that can be computed from EnergyPlus or metered from real devices. In fact, summing up the delivered heating/cooling to all zones, 's will not necessarily give the total electrical energy consumption of the building. Therefore, in this paper, we will use a combined statistical and EnergyPlus approach to estimate the actual energy consumed by the building. The approach will be discussed in Section IV.
We estimate the building total power as a function of and with the latter one usually having sufficient historical data. Let us denote by the total power at time ; then the relationship between aforementioned variables can be shown by , where is a vector of 's and is a vector of ones. Next, we will present how to estimate using a generalized form of Cochrane-Orcutt estimation technique.

IV. ENERGY FORECAST MODEL
The relationship between , total power at time and can be modeled through a simple linear regression as follows: (13) where is the error term at time and 's are linear model parameters. If the assumption of linearity were met, then would typically be assumed independent and the model parameters, 's, would be estimated using Least Squares Estimation (LSE) technique. However, the actual relationship between total power consumption at time with effective cooling power and external temperature may follow an unknown nonlinear model. In addition, there are more variables such as occupancy and cooling fans power that can affect the total power consumption. These effects cannot be explained through the linear structure of (13), and as a result, they emerge into the error terms. In this situation, the assumption of independency is no longer met and the ordinary LSE technique cannot be applied [22]. To avoid this problem, we employ the Cochran-Orcutt technique by rewriting (13) as follows: (14) is an independent and identical white noise and is a function of past error terms representing the structure of an autocorrelation. Following [24], we obtain transformed variables: (15) where and are autoregressive operators with orders of and that are applied to both the external temperature and the vector of zonal effective power to find the building total energy demand. Hence, (14) is replaced by (16) This is a typical linear regression with independent error terms and parameters that can be estimated iteratively as proposed in [24].

V. OPTIMAL CONTROL STRATEGY
In this section, we propose an optimal control strategy that includes developing a mathematical programming formulation that is solved dynamically over time. The cooling/heating model presented in (6) is a dynamic model that describes how the state variables, 's, evolve over time by starting from an initial condition and by manipulating control variables, 's. At time , the actual internal temperature, is unknown and is specified by replacing it with any arbitrary set point value, i.e., . Then, can be calculated accordingly using (7). is then fed into (13) to calculate the corresponding building total energy use. This is repeated for the next 24 h and for all combinations of set point values between and and all zones to find the optimal combination of set points that minimize total energy use and total deviation from the thermal comfort.
The dynamic model requires a dynamic programming scheme to find the optimal control variables in such a way that the objective function is optimized over a specific time horizon . In this study, we set , so that the control scheme can provide the optimal control variables for any 24 look-ahead periods (hours in the current case). At any given hour, the optimization procedure is repeated for the next 24 h by updating the state variables 's, and external temperature values. At each time step, the optimal control scheme solves the following multi-objective problem: subject to the following constraints: where and are the thermal comfort upper and lower bounds for the th zone internal temperature at time . is the temperature violation below the lower bound and is the temperature violation above the upper bound. is the internal temperature at time 0 and is the length of the time period that is set equal to 1 h. The thermal comfort constraint is imposed to each single zone based on the current zone temperature. There are two objective functions: (1) is the total cost of energy, which includes Total Usage Cost (cost per kWh) and Daily As-Used Demand Charges. The latter cost is determined for each weekday in the billing period and applied to the daily peak demand during each time period. The Monthly As-Used Demand Charges for the billing period are equal to the sum As-Used Daily Demand charges for the time periods [28].
Here, is unit price of electricity at time and is the total penalty for exceeding the thermal comfort bounds. This function penalizes any deviation from the predefined comfort bounds at any given time. is the penalty applied on peak energy demand, is a set representing the on-peak period, and is the penalty that is applied to any violation from the comfort bounds at time . The latter parameter indicates different discomfort costs for different hours of a day. Since and do not match in units and scale, it is not possible to integrate them into a single objective function by simply adding their values. Hence, we build a Multiobjective Mathematical Programming (MMP) structure using a weighted metric method. This method requires less restrictive assumptions (see, e.g., [25]) compared to other methods. This method scales and into a single objective function that is in an metric form as follows: (17) Both and in and are replaced with (17). In this equation, and are the non-negative weights for and , respectively, such that . and represent the relative importance of the objective functions and are determined by the decision maker.
indicates the type of metric we use in our problem. For example, is equivalent to solving the weighted sum of deviations from ideal values, whereas, means minimizing the weighted Euclidean distance of any point in the objective space from its ideal point. It is proven that if there exist bounded solutions for and , then for any combinations of 's and values, there is one or more Pareto solutions [25], [26]. A Pareto solution has the property that, for any point in the Pareto set, there does not exist another point that results in better performance for both objectives simultaneously [26]. and are the minimum and maximum possible values of the th objective function. It is easy to find the minimum and maximum possible values for the second objective function. The minimum penalty for exceeding the thermal comfort bounds is 0, when the internal temperature is within thermal comfort or equivalently . The maximum penalty for exceeding thermal comfort, , can be calculated by multiplying the maximum penalty, and the maximum deviation form thermal comfort as follows: (18) Similarly, the minimum value of the first objective function, is set equal to zero. This is because the first objective function is the total energy cost, and when there is no demand, the cost is zero.
can be found using EnergyPlus by calculating the maximum energy consumed by HVAC. That is (19) where is the maximum load that can be generated by HVAC at any time and is the highest external temperature. Note that (19) calculates the maximum cost for the maximum energy consumed by HVAC for the case of the highest daily temperature. This gives an upper limit for total energy cost at any time.
Both components of (17) are normalized between 0 and 1 and as a result the metric function varies between 0 and 1. One disadvantage of the weighted metric method is that its performance highly depends upon the values of and . For example, the maximum HVAC capacity and the highest external temperature may result in an extremely large value of , that is, the first component of (19) becomes much smaller than the second component (thermal comfort). Subsequently, the thermal comfort influences the optimization results more than is desired. This undesirable effect can significantly be reduced by incorporating domain specific knowledge on . Replacing and with in (17), the model now comprises an ordinary mathematical program with a single objective function that can be solved using Dynamic Programming. At time , the model is solved and the optimal control values for the next 24-h are obtained. Then, in the next hour, once some new data from each of the zones' internal temperature are received, we update all of the state variables and solve the model with the updated forecast values. This process is repeated every hour and once new information is received, the optimization process will repeat for the next 24 h.

VI. NUMERICAL EXAMPLE
In this section, we evaluate the performance of our proposed control strategy using an illustrative large size office building. This is a commercial reference-building developed by the U.S. Department of Energy (further details can be found in [27]). In Offline Phase (see Fig. 1), we select 12 main zones of the building and run the simulation using 2009 Phoenix weather data. We collect one-month of hourly data (July 2009) and extract 744 hourly basis simulated output data for the following variables: zonal internal temperatures, building external temperature, zonal cooling power rate values, and the building total energy demand. As mentioned in Section II, in Online Phase, the proposed model can be updated using actual data. However, in this study, since we do not access to such information, we again use simulation model to generate required data over time. We mimic the Online Phase (see Fig. 1), by using one-month of simulation data for August 2009. For each hour, the simulation model is run once, the model parameters and optimization model are solved, and optimal set point values are updated. Phase 2 outputs are: the total energy cost including usage cost, daily as-used demand charges and thermal discomfort penalty, and the set of optimal set point values for each zone at time to time . The estimated parameters of the proposed cooling/heating model are presented in Table I. Without loss of generality, we assume that and are constant over time or equivalently and for . It implies that during simulation the overall weather pattern does not radically change, and there is no latent variables (such as heat gain/lost by lowquality insulation or other unknown parameters). This assumption may work for a short time; however, for a long-term optimal control strategy one should estimate the parameters adaptively (e.g., as a function of measured or forecast weather conditions). As shown in Table I, 's are all showing negative correlations between the next hour internal temperature and zonal effective power. This is because the simulated dataset is run for a summer month and the effective power is used for cooling of the building. Similarly, 's show negative correlation between current internal temperature and external temperature. Fig. 2 displays a radar chart for the correlations between the actual internal temperature values and their corresponding -step-ahead forecast values for each zone. In this figure, the th radius represents the th zone and the th irregular polygon represents the correlation between the -step-ahead forecasted internal temperature and its corresponding actual temperature. It shows that the correlations between actual internal temperature values and their corresponding one-step-ahead forecast  values vary between 0.75 and 0.85. This also shows that the proposed cooling/heating model can appropriately provide the one-step-ahead forecast for the zonal internal temperature. Furthermore, it is observed that the performance of the proposed heating/cooling model decreases as the lag increases. For instance, the correlation between the actual and forecasted internal temperature for the third lag varies around 0.5, which is expected. This also can be seen in Fig. 3. In this figure, we plot the first zone's actual and forecasted internal temperature values for the first and fifth lags. The one-step-ahead forecast values are quite close to the actual data while the five-step-ahead forecast values are relatively far from the actual internal temperature. As previously discussed, this is mainly because we assume that the parameter values do not vary during every 24 h or equivalently and for . It means that the structure of balance equations does not change as time passes. The performance of the proposed forecast values can be improved by applying a time-dependent parameter estimation method.
To evaluate the performance of the energy forecast model we use simulated data. Table II presents the estimated parameters of  TABLE II  ESTIMATED PARAMETERS OF ENERGY DEMAND FORECAST MODEL   TABLE III  for the 24th lag are 95.9% and 95.7%, respectively. This indicates that the proposed model can provide high quality 24-h-step-ahead forecast values for total energy demand, which are very close to the actual energy demand. Fig. 4 presents a comparison between the actual and forecasted total energy demand for the 24th lag. This forecast model can well capture the majority of the variations. Thus, the proposed forecast model is deemed a reasonable replacement for the simulation, particularly well-suited for optimization purposes.
We now investiage the performance of the proposed optimal control strategy for the above illustrative model. Table III presents the input parameters of the proposed mathematical model. Note that , the thermal discomfort penalty is a normalized multiplier, defined by the decision maker (e.g., building energy manager), to weigh the importance of keeping the internal temprature within the thermal comfort at time .
is normalized based upon the number of people working in each zone such that the total thermal discomfort penalty over 24 h is equal to 1.0.
We assume that the internal temperature values can vary between 62 F and 76 F with an increment rate of 0.5 units. However, and in are fixed and set equal to 68 F and 72 F, respectively. Therefore, any internal temperature value less than 68 F or greater than 72 F is penalized by applying a thermal discomfort cost. Besides, and , the importance factors associated with objective functions and are set equal to 0.7 and 0.3, indicating that the energy cost is relatively more important than the thermal discomfort penalty. Without loss of generality, we assume , which turns the problem into a weighted sum of normalized deviations. Fig. 5 illustrates the effect of the proposed weighted metric method in combining both objective functions. The figure plots the minimum cost-to-go values (minimum total cost from step 1 to step ) in dynamic programming. The axis shows the feasible set point values ranging between 62 F to 78 F and the axis shows the minimum values of the combined function, , over the reduced horizon from to when . For the first 6 h, as set point values increase, the minimum cost-to-go value decreases. This is because we set the thermal discomfort penalty equal to zero for the first 6 h (see Table III). This implies that for the first 6 h, the only active objective function is . In this case, higher temperature values result in lower demands without increasing the thermal discomfort objective function. Fig. 6 presents values of optimal set point, energy demand, and operating costs versus time for a period of 72 h. Fig. 6(a) for instance, depicts the optimal set point values and their corresponding values of internal temperature.
As shown in Fig. 6(a), the optimal set point values are very low just before the on-peak period begins. This means that the proposed control scheme suggests precooling the building before the price goes up and before the daily as-used demand charge becomes effective. Then during the first hours of the on-peak period, the cooling stored in the zone thermal masses is discharged. This helps the building save the most possible energy during peak time.
It is worth noting that in an ideal HVAC system, the optimal set points and internal temperature are same values. However, in the real world, many unknown or uncontrollable factors can significantly cause the actual internal temperature to be deviated from the set point values. This is particularly obvious when the proposed control scheme select a very low set point at time after a higher set point value at time [e.g., see in Fig. 6(a)]. In this situation, HVAC system cannot reach the optimal set point value in such a short time frame. In fact, this is one reason that the model is updated when new information from internal temperature is received. Fig. 6(b) illustrates that total building energy demand oscillates according to the zone effective power and set point values. Fig. 6(c) includes separate profiles for total use cost per kWh, the daily as-used demand charge cost and total energy cost. Note that in this example, it is assumed that the daily as-used demand charge is applied once in a day to the maximum energy consumed in peak hours between 1:00 p.m. and 6:00 p.m. Therefore, the precooling time is immediately before 1:00 p.m. and is applied to lower the maximum energy demand after 1:00 and as a result reduce the daily as-used demand charge cost. Fig. 7 presents further details about how the two different objective functions are combined to shape the single measure that is optimized. Fig. 7(a) is the total energy cost including energy costs and daily as-used demand charges. Fig. 7(b) depicts the temperature violation from the thermal comfort bounds over time. As shown in this figure, the maximum and minimum values for each of these functions are different. Although the violation from thermal comfort is slightly higher in off-peak periods, when less people are in the building, there are still large violations in on-peak periods. This is because, in on-peak periods, the energy demand is very high and the proposed model  attempts to reduce the costs by allowing some violations above or below the thermal comfort bounds. A normalized version of this figure is presented in Fig. 8. These plots present the normalized profiles associated with total energy cost, , as well as the total thermal comfort violation penalty . The maximum values of these two functions are slightly different. The maximum values of are in peak hours while the maximum values of occur immediately before peak hours. This configuration leads to the minimum combined metric profile that is . We also compare our proposed control scheme with a number of simple alternative control schemes, as shown in Fig. 9. These configurations are set to cover different scenarios and provide good insights about the proposed control scheme. The "Constant Controller 1" and the "Highest Value" control schemes have fixed set point values that do not change during the 24 h period shown. The "Highest Value" controller is set equal to be 76 F and as a result, returns the lowest total energy cost but the maximum thermal discomfort cost for this particular day that requires cooling. The dynamic controllers vary over time using different patterns. For example, "Dynamic Controller 2" allows HVAC to provide more cooling before peak hours and then during peak hours it increases the set point values. We use different controllers to understand the characteristics of the proposed control scheme. Table IV depicts the results of running different control schemes for the medium office building model. It shows that the proposed scheme is superior to all of the other alternatives in the total combined metric value, i.e., . This means that the proposed control scheme is able to find the best compromise between energy costs and thermal comfort. The energy cost is minimal for the "Highest Value" controller. This is obvious, because it provides the minimum possible cooling load and keeps the internal temperature around 76 F. However, the normalized thermal discomfort penalty , for such controller is 0.81, which is significantly greater than the other controllers. Such an extreme controller provides internal temperature values that are often beyond the upper thermal comfort bound and are impractical in the real world. Note that the maximum thermal discomfort penalty, , is designed to be 1, where the internal temperature would never go below the maximum thermal discomfort value. To calculate average daily violation from Thermal Comfort (F), we calculate the deviation from thermal comfort for each hour at each zone, then we take the average of these deviations over 24 h for all zones. The normalized energy cost, for the constant controller (Set F) is 0.27, which is the highest energy cost among the alternatives. The minimum normalized energy cost is associated with the proposed control scheme (and Dynamic Controller 2) and is equal to 0.20. This implies that the proposed control scheme offers the set point values with minimum cost while maintaining the thermal comfort in an appropriate level, thus is minimized. Although the energy costs for the various control schemes are clearly different, their normalized values of energy cost are very close to each other (from 0.20 to 0.26). This is because the maximum value of energy cost is very large. If is large, then the first component of (17)-the normalized energy cost-becomes very small. To tackle this problem, one can propose better choices for or apply other multi-objective techniques to combine the two objective functions. Table IV reveals that the proposed framework can adequately be employed to minimize building energy costs and at the same time keep the internal temperature in the range of thermal comfort bounds.

VII. SUMMARY AND CONCLUSION
The above results can be summarized as follows: (i) The proposed heating/cooling model can appropriately forecast the internal temperature values for the next few hours, but the methodology needs to be improved for longer lags. (ii) The energy forecast model works quite well for total energy demand. (iii) the proposed control scheme provides the set point values that minimize both total energy cost and thermal discomfort penalty at the same time. To improve the performance of the proposed heating/cooling model, one can develop an adaptive approach-Bayes estimation for instance or considering longer time history in (6)-to estimate the model parameters. Currently, the forecast errors are adjusted once new measured data is received. However, the heating/cooling model does not utilize any new information to update the parameters and forecast values, i.e., it is not currently an adaptive model. Also, additional zonal variables such as occupancy, heat transfer due to inter-correlation between zones, and zone surface phenomena can be added into the proposed heating/cooling model to improve the quality of forecast values in farther lags. Another potential improvement for the proposed control strategy is to employ more efficient methods for the MMP problem. The Utility Function approach, for instance, can be used instead of the weighted metric method to determine the lowest and highest desirable values of the objective functions.