Building energy consumption on-line forecasting using physics based system identification

Energy and Buildings 82 (2014) 1–12

Contents lists available at ScienceDirect

Energy and Buildings journal homepage: www.elsevier.com/locate/enbuild

Building energy consumption on-line forecasting using physics based system identification Xiwang Li ∗ , Jin Wen Department of Civil, Architectural and Environmental Engineering, Drexel University, 3141 Chestnut Street, Curtis 251, Philadelphia, PA 19104, USA

a r t i c l e

i n f o

Article history: Received 14 March 2014 Received in revised form 8 July 2014 Accepted 9 July 2014 Available online 18 July 2014 Keywords: Building energy modeling System identification Building control and operation Building energy efficiency

a b s t r a c t Model based control has become a promising solution for building operation optimization and energy saving. Accuracy and computationally efficiency are two of the most important requirements for building energy models. Existing studies in this area have mostly been focusing on reducing computation burden using simplified physics based modeling approach. However, creating even the simplified physics based model is often challenging and time consuming. Pure date-driven statistical models have also been adopted in a lot of studies. Such models, unfortunately, often require long training period and are bounded to building operating conditions that they are trained for. Therefore, this study proposes a novel methodology to develop building energy estimation models for on-line building control and optimization using a system identification approach. Frequency domain spectral density analysis is implemented in this on-line modeling approach to capture the dynamics of building energy system and forecast the energy consumption with more than 90% accuracy and less than 2 min computational speed. A systematic analysis of system structure, system order and system excitation selection are also demonstrated. The forecasting results from this proposed model are validated against detailed physics based simulation results using a mid-size commercial building EnergyPlus model. © 2014 Elsevier B.V. All rights reserved.

1. Introduction Buildings roughly account for 40% primary energy use in U.S. It costs a total utility bill of more than $400 billion in 2012. Around 30% of the energy used in building is consumed by heating, ventilating and air conditioning (HVAC) [1]. Study has shown that non-optimal control of building energy system could cause the malfunction of equipment or performance degradation from 15% to 30% in commercial buildings [2]. Moreover, it is estimated by the National Energy Technology Laboratory that more than onefourth of the 713 GW of U.S. electricity demand in 2010 could be dispatchable if only buildings could respond to that dispatch through advanced building energy control and operation strategies and smart grid infrastructure [3]. Therefore building control and operation is significant economically and environmentally. In strategies used to improve building control and operation, high fidelity building energy model is the most critical component. Existing building energy models can be categorized into white box (physics-based) models, black box (data-driven) models and grey box (hybrid) models. One of the most comprehensive white

∗ Corresponding author. Tel.: +1 215 895 6941; fax: +1 215 895 1364. E-mail addresses: [email protected], [email protected] (X. Li). http://dx.doi.org/10.1016/j.enbuild.2014.07.021 0378-7788/© 2014 Elsevier B.V. All rights reserved.

box models is EnergyPlus, which is a whole building energy simulation program that engineers, architects, and researchers use to model energy and water use in buildings [4]. A Building Controls Virtual Test Bed (BCVTB) was developed by Wetter and Haves to link the building models with control systems [5]. BCVTB is a middleware tool that allows to data sharing among different simulation programs, such as EnergyPlus, Matlab and Modelica, for distributed simulation. Therefore, through this test bed different user defined building control and optimization strategies can be applied into different building simulation models. Ma et al. proposed and demonstrated an economic MPC technique to reduce energy and demand cost [6]. Using BCVTB as middleware, realtime data exchange between EnergyPlus and a Matlab controller was realized. An economic objective function to minimize daily electricity costs was then developed in MPC and applied in EnergyPlus model through the middleware. About 25.3% energy saving and 28.5% cost saving were achieved by this MPC in a single story commercial building located in Chicago, Illinois. Corbin et al. [7] utilized a Matlab-EnergyPlus MPC environment and incorporated it with a particle swarm optimizer to predict optimal building control strategies. Every time step, the schedules and setpoints will be wrote into EnergyPlus input IDF file and then EnergyPlus output results will be evaluated within the MATLAB optimal control module, based upon the objective cost function, then the

2

X. Li, J. Wen / Energy and Buildings 82 (2014) 1–12

updated temperature setpoints will be sent to EnergyPlus again. This procedure will be repeated until some convergence criterions are satisfied. This on-line optimization environment was also applied in a DOE benchmark building EnergyPlus models [8]. The results showed 5% cost saving just by optimizing the hourly cooling setpoint in large office building model, and 54% energy saving by determining hourly supply water temperature in small office building model. Another important control optimization software environment, Genopt, was also developed by Wetter in 2000 [9], which can iteratively execute any simulation program based on plain text input/output files until an optimal solution is found. Genopt was used by Coffey et al. [10] to incorporate a modified genetic algorithm MPC with an EnergyPlus model to study the temperature control optimization in office buildings and its effect building energy demand. Even though these elaborate simulation tools are very effective and accurate, they require detailed information and parameters of buildings, energy system and outside weather conditions. Identifying these parameters, however, is very time consuming and need expert work. What is even more challenging, the simulation speed is relative low and not suitable to be used in on-line MPC. Therefore, the two major researching objectives of recent studies on building MOC are increasing simulation speed and maintaining simulation accuracy. Cole et al. [11] used OpenStudio to reduce the EnergyPlus model by perturbing the system and fitting the results into a reduced-order linear model The reduced model is able to simulate an one-day simulation within 1 s and the discrepancy between this reduce-order linear model and EnergyPlus model is less than 2.3% under an optimal temperature setpoints setting situation. Black box model is also known as purely data driven model. Statistical models are simply applied to capture the correlation between building energy consumption and several building operation variables. This method needs the on-site measurements over a certain period of time to train the statistical models which can calculate the accurate predictions under different conditions. Autoregressive with exogenous (ARX) model was implemented to predict the 1 h ahead building load by Yun et al. [12]. This predictive model is applied on several different DOE benchmark buildings [8] to choose the building control strategies. Artificial neural network (ANN) is another popular method in building energy prediction for building operation purpose. Adaptive ANNs for on-line building energy prediction was investigated by Yang et al. [13]. These models are capable of adapting themselves to unexpected pattern changes in the incoming data, and therefore can be used for the realtime on-line building energy prediction, which is the fundamental for building optimal control. Chen et al. [14] developed a day-based wavelet ANN method for next day load forecasting. They selected a historical day with the same weekday index and similar weather condition as next day. They similar day load was then decomposed into multiple levels by using wavelet decomposition. These decomposed components were fed into different ANNs to predict the next day load at each component. After all, all these prediction results will be combined into an overall forecasting result. Black box models are easy to build and computationally efficient, however, such models often require long training period and are bounded to building operating conditions that they are trained for which sometimes can cause huge forecasting error when training data does not cover all the forecasting range, and moreover, black box always sacrifice physical insights to obtain high accuracy and calculation speed. Resistance and Capacitance (RC) network model is the most common grey model for building energy estimation, which requires less training data and has less parameter to determine. RC model usually uses state space model to estimate the building heating and cooling load [15], and control building temperature [16,17]. Different optimization and searching algorithms have been utilized in determining the resistances and capacitances [18,19]. Ma et al. [20]

reduced the order of RC model to increase the calculation speed, sing balance realization method. It still took 28 min for this reduced model to run a one-week simulation using a reduced RC model, which is still not suit for on-line building operation optimization. Different to pervious works using single RC models, an integrated 3R2C and EKF (Extended Kalman Filter) model was developed to estimate the building energy consumption in Ref. [21]. In this work, an EKF was used to estimate the state vector X using real sensor measurement data. The estimated load matched the EnergyPlus results within 10% at 93% of the time. It is common knowledge that determining the parameter of RC model is computational demanded and the model structures and parameters are unique from building to building. It is impossible to utilize a RC model for one building to another building. System identification is the process of developing or improving a mathematical representation of a physical system using experimental data [22]. System identification techniques have started been applied in building energy model area to obtain better and more accurate estimation of building performance. Privara et al. [23] proposed an approach combining a EnergyPlus model and a subspace system identification model to forecast the building performance. In this study, a six-floor office building was modeled by a Matlab toolbox called N4SID to estimate the building zone temperature. However, directly using Matlab toolbox cannot guarantee the forecasting accuracy every time. The accuracy of Matlab toolbox is depended on the data properties. It is very hard to guarantee the toolbox suitable to every building under different operation strategies. On the other hand, there is no systematic discussion of model structure selection and parameter determining presented in this paper. Agbi et al. [24] studied system structural identifiability and model parameters identifiability for building energy state space model. Unfortunately, this paper focuses on the discussion of identifiability of the system upon available data without any active building excitation. It did not improve the system identification model developing speed or improve the model forecasting accuracy. And it still uses the RC network model, so the method demonstrated in this paper has the limitation of all the RC model has. For example, the parameters identification for those Rs and Cs is very difficult and time consuming. In this study, an EnergyPlus model of a mid-size commercial building was used to provide building operation data in lieu of a real building. Then a systematic discussion system identification model development and system excitation base on those operation data is presented in this paper. The combination of active system excitation and a frequency domain spectral density analysis are applied to develop building energy forecasting models that are genetic (can be easily extended to other building types/operating conditions), accurate, and computationally cost-effective. The system identification model development will be discussed in Section 2, followed by the development of system identification based building energy on-line forecasting model in Section 3. The forecasting results will then be analyzed in Section 4. After this system identification approach has been developed and validated in a small commercial building, which is then applied in a medium building.

2. Building energy model development In this study, the on-line building energy forecasting model is developed based a mid-size reference commercial EnergyPlus simulation model, provided by U.S. Department of Energy (DOE) [25]. This reference building EnergyPlus model has been validated and used for optimize designs, operation and advanced controls in numbers of studies [26–29] The procedure of model structure determining, input and output selection, system exciting and model training and validation will be introduced in detail in this section.

X. Li, J. Wen / Energy and Buildings 82 (2014) 1–12

3

Table 2 Variables of system identification model. Parameter zone model

Fig. 1. Small commercial building view.

2.1. Building description This single story office building (Fig. 1) has five zones, and the total floor area is 510 m2 . The window-to-wall ratio of this building’s facades is approximately 21.2%, and the windows are equally distributed. The U-factor of these single pane windows is 3.4 W/m2 K and the solar heat gain factor is 0.36. The solar absorptivity, transmissivity and reflectivity are 0.06, 0.697 and 0.243, respectively. The location of this building selected for this study is Philadelphia, PA, USA. This building is partitioned into five different air conditioning zones, and an unconditioned attic zone. The five conditioned zone are four perimeter zones and one core zone. The roof insulation has an R-value of 15. The roof is covered in an asphalt membrane, with a solar absorptivity value of 0.9. The exterior wall has the following construction (from outside layer to inside layer): • • • •

2.54 cm of stucco. 20.32 cm concrete masonry units. R-6 continuous insulation. 1.27 cm gypsum wallboard.

Since at the purpose of this study is to develop a building energy forecasting model using system identification approach. For the simplification of the mechanical system, and mainly focusing on the building thermal response dynamics, the mechanical system in this building is constant-air-volume (CAV) air handling units (AHUs). Direct expansion (DX) coil is used in this building, which is connected to an electrical cooling source with coefficient of performance (COP) of 3. Heating is provided by gas boiler with efficiency of 0.95. The baseline model internal load inputs are summarized in Table 1. All these assumptions are from Deru et al. [25]. These model assumptions will be updated in the system identification process in order to enrich the operation data. 2.2. System identification model selection As stated above, the objective of this study is to develop an on-line building energy model using system identification method. Model structure plays the most important role in model forecasting accuracy. Privara et al. summarized and compared different system identification approaches for building energy modeling and control [30]. There are two different categories of system identification models, which are time domain models and frequency domain models. There is large number of modeling and Table 1 Building model internal and external gains. Variable

Value

Occupant density Ventilation requirement Lighting power density Interior small plug loads Envelope infiltration rate

0.1 Person/m2 2.36e−3 m3 /(people m2 ) + 3.05e−4 m3 /(s m2 ) 18 W/m2 10 W/m2 9.4e−4 m3 /(s m2 )

Variable

Variable name

Type

Qload Ezone Tout Tzone Tadj Rin Tsol-air Qvent Qdir Qdif Qfan

Zone heating/cooling load (W) Zone heating/cooling energy (J) Outdoor air temperature (C) Zone temperature (C) Adjacent zone temperature (C) Equipment and occupancy schedule ratio (–) Solar air temperature (C) Ventilation rate (m3 /s) Direct solar radiation (W/m2 ) Diffuse solar radiation (W/m2 ) Supply fan heat gain (J)

Output Output Input Input Input Input Input Input Input Input Input

identification approaches developed over the recent years, but few of them are suitable for building energy modeling, neither a universal model that can work for different building types and operation schemes. Building energy systems, especially the HVAC systems, are very complicated nonlinear dynamic systems, which is hard to model and forecast. In order to develop a relative simple system model, frequency response function approach is applied in this study due to its excellent performance in handling system nonlinearity. Fundamentally a frequency response function is a mathematical representation of the relationship between the input and the output of a system in frequency domain, which is a special case of the transfer function, but it can simplify the time domain transfer function and maintain the useful information (Eq. (1)). H(jω) =

Syu (jω) Y (jω) = U(jω) Suu (jω)

(1)

where Y(jω) the Fourier is transform of system output y(t), and U(jω) is the Fourier transform of system input U(t). However, better results can be obtained in practice by computing the frequency response function (Syu ) as the ratio of cross-spectrum between input and output to the power spectrum of the input (Suu ). Then by applying the Inverse Fourier Transform, the Impulse Response Functions (IRF) per measurement channel are obtained. In this way, the operational data across the full spectrum of the building operation range will be analyzed accurately. The detailed information about the spectral density model implementation in this study will be presented in Section 3. In this project, two different energy models will be developed for the core zone and perimeter zone of this building separately, which is because the outside weather condition directly affects the behavior of perimeter zone but not to core zones. Therefore, the perimeter zone model needs to include outdoor weather information, such as temperature, solar radiation, etc. in the input. 2.3. Model input and output determining Besides the model structure, model input and output selection is also crucial to the accuracy of system identification model. The input and output selection should be based on the physical relationship and the data availability. From the objective of this system identification model, two different models have been developed. The outputs of the model are relative easy to determine according to our modeling objective. They are building cooling and heating energy consumption or building heating and cooling load. The inputs and outputs of perimeter zone system identification model are tableted in Table 2. The inputs and outputs of core zone model are very similar to those of perimeter zone model, but exclude Tsol-air , Qdir and Qdiff in system inputs, because there is no direct solar radiation or transmission in core zone. Solar air

4

X. Li, J. Wen / Energy and Buildings 82 (2014) 1–12

temperature (Tsol-air ) is a variable used to determine the total heat gain through opaque exterior surfaces to calculate cooling load of a building. It is not an output from EnergyPlus model, nor a measurement from weather station, and then a specific calculation model is created, using the following equation: Tsol−air = To +

˛I − Qir ho

(2)

where ˛ is absorptivity of each opaque wall; I is global solar irradiance (W/m2 ); and Qir is extra infrared radiation due to difference between the external air temperature and the apparent sky temperature (W/m2 ). Direct solar radiation and diffuse solar radiation are used to estimate the building heat gain due to the solar transmission through windows. They can be either obtained from weather forecasting information or calculated from global solar irradiance. All these input variables are categorized into two groups: unexcited and excited inputs. The zone temperature, adjacent zone temperature and equipment and occupancy schedule ratio are excited inputs, which will be excited during model training period, and the other variables are unexcited inputs. The zone temperature is excited by changing the zone temperature setpoints, and the equipment and occupancy schedule ratio is excited by updating their on/off schedule ratio in EnergyPlus model. The details about the excitation signal generation will be introduced in Section 2.4. 2.4. System excitation signal generation In order to develop the spectrum density model and improve the accuracy of system identification model, exciting signals were generated and injected into the system during the training period. These exciting signals include zone temperature setpoints, internal equipment and occupancy schedules. Such signals are needed to satisfy key theoretical assumptions on reliable statistical identification – persistent exciting signals [31]. Three basic “plant-friendly” excitation signal constraints were presented by Braun et al. [32]: • Keeping minimum deviations in product quality; • Implementing a signal of short duration to minimize the amount of off-spec product; • Keeping move sizes small to satisfy actuator constraints and minimize “wear and tear” on process equipment. Therefore, certain constraints added into our system excitation process, for example, the boundary of temperature setpoints, the minimum temperature setpoint or equipment schedule updating time span, and so on. Different system excitation strategies, such as Pseudo-Random Binary signal, Pseudo-Random Sequences, Multisine signal, have been discussed and applied in on-linear process systems in different areas [32–34]. However, there are few publications on building energy system excitation found in existing literature. Pseudo-random Binary Sequence (PRBS) excitation signals (Eq. (3)) for building temperature setpoints was generated and applied in [6].

⎧ ⎨ 21, excited zone, PRBS = 0 Tsp,i (k)

⎩

25, excited zone, PRBS = 1

(3)

25, non-excited zones

2.4.1. Excitation signal generation function The objective of this on-line model is for building operation and control, therefore, the excitation signals do not need to cover the entire frequency domain, but the building operation range. To guarantee those signals cover the frequency domain around the building

Fig. 2. Excitation signal generation procedure.

common operation range, sum of sinusoids (SINE) model was used to generate the exciting signal (Eq. (4)). U+1 = U + A



2˛j sin(ω tT + ϕ )

(4)

where U+1 is the excitation signal; A is a magnitude scale parameter from 0 to 1; ω is periodic frequency parameter from 0 to 2; T is the sampling time, ϕ is the phase lag parameter from 0 to 2, and ˛j is the Fourier coefficients. These parameters will be updated every time step. Another benefit of using sum of sinusoids input signals is that they enable the user to directly specify the shape and character of the power spectrum. The guidelines of excitation signal generation function parameter determining from Rivera et al. [35] is applied to ensure the signals contains necessary frequency information: 1 H ˇs dom

≤ω≤

ϕ = 2

i 

˛s L dom

j˛j

(5)

(6)

j=1 H dom

L and dom correspond to the high and low estimates of where the dominant time constant of the system (denote the slowest and the fastest systems time constants). ˛s and ˇs are user-decisions on high and low frequency content based on identification requirement. Typically, ˛s is 2 and ˇs is 3, corresponding to 95% of settling time [35]. In additional to function magnitude scale parameter, periodic frequency parameter and phase lag parameter, harmonic suppression is another important parameter which can decompose the output signal to obtain a more accurate estimate of the linear and nonlinear components.

2.4.2. Excitation function parameter specification and data generation The procedure of excitation signals design is summarized in Fig. 2. The response time constant of a dynamic system is a measure of how quickly the system responds. It is usually measured by experiments. For example, if the impulse response of a dynamic system can be expressed as: x(t) =

˛ T

e−t/T

(7)

where T is the response time constant, ˛ is a state parameter. When x(t) = 0.95, t = T.95 . Therefore, with considering the building T dynamics of the target building in this project, dom = 360 min, and L = 30 min; ˛s and ˇs are 2 and 3, respectively. For the temperdom ature setpoint excitation signals in this project, Tmax = 32 ◦ C(90 ◦ F) and Tmin =10 ◦ C (50 ◦ F); while for the schedule ratio excitation signals, Rmax = 1 and Rmin = 0.

X. Li, J. Wen / Energy and Buildings 82 (2014) 1–12

5

Table 3 Excitation frequency and sampling length testing summary. R2

Excitation frequency 15 min

Training period Forecasting period

Sampling length 3h

Excitation frequency 30 min

93.65% 95.30%

Once all the time independent parameters are determined, the excitation function will generate the signals every time step, which will then be applied in building model automatically. 2.4.3. Excitation signal injection frequency According to Braun’s guidelines, the excitation injection frequency cannot be too high. The injection time step should be larger than the system response time. The building response time is related to building thermal mass, which can only be found through experiment test. The building studied in this project is a small light building. It system response time is relatively short. A parametric experiment test has been conducted to find out the best injection frequency. However, once the system excitation frequency changed, the building dynamics would change, and the power density model sapling length should also be updated accordingly. Table 3 presents the excitation frequency and sampling length testing results. Finally, 30 min with 6 h and 30 min with 4 h have been chosen for excitation frequency and sampling length for core zone model and perimeter zone model, respectively. One case temperature setpoints and equipment schedule excitements used in this study are shown in Fig. 3, where all the parameters are updating every time step and the excitation signals are injected into the building model every 30 min. 2.5. Training and validation data generation The system excitation signals discussed in Section 2.4 was modeled and generated in Matlab. In order to apply these excitements into EnergyPlus model, BCVTB is used to exchange data between Matlab and EnergyPlus. Here BCVTB plays a master role in data exchange between Matlab and EnergyPlus through runtime coupling, as shown in Fig. 4. During the entire study, typical meteorological year (TMY) weather data for Philadelphia is used. During the training and validation period, excited and unexcited building control signals will be sent to EnergyPlus model following the procedure in Fig. 2, respectively. Simulation results and control signals will be sent back and stored in Matlab for system identification model training and validation, as discussed in Section 2.3.

Sampling length 4h

93.70% 95.30%

Excitation frequency 60 min

Sampling length 6h

92.21% 91.01%

The connection between Matlab and EnergyPlus through BCVTB is illustrated in Fig. 5, where the Matlab simulator is the excited or unexcited control signal generator, and EnergyPlus simulator is the building model discussed in Section 2. The length of training time is changeable according to the forecasting accuracy requirement. 3. Building energy on-line estimation model 3.1. Building energy system identification model development The system identification model structure selection has been discussed in Section 2.2. Frequency response function model has been chosen for this project. Fig. 6 shows the system identification model development from building operation data using spectral density model for frequency response function. In this figure, U is training inputs, h is a reference signal to analyze the input data, Y is training outputs data, PSD is power spectral density model for inputs data, CPSD is cross power spectral density model for input and output, Suu is the result of PSD, Syu is the result of CPSD, G is the transfer function from input to output, and yˆ is the output estimation. Suu and Syu are estimation the correlation between input and output. G(z) is the transfer function in frequency domain, which can be transferred to time domain transfer function G(t) using inverse Fourier function transformation. As discussed before, in order to calculate the frequency response function, power spectral density model is applied. Power spectral density describes how the power of a signal or time series is distributed over the frequency spectrum, which is property of the system signal and very useful in frequency domain system identification (Eqs. (8) and (9)). 1 −j 2k l Ruu ()e l

(8)

1 −j 2k l Ryu ()e l

(9)

l−1

Suu =

=1 l−1

Syu =

=1

where Ruu is the auto-correlation between the inputs and Ryu is the cross-correlation between input and output [22], which are calculated in Eqs. (10) and (11). u is the input vector, y is the output

Fig. 3. Temperature setpoint and equipment schedule excitement.

Fig. 4. Building operation data for on-line model training and validation.

6

X. Li, J. Wen / Energy and Buildings 82 (2014) 1–12

Fig. 5. Matlab-BCVTB-EnergyPlus model.

vector, and l is the length of the sampling data. l is a very important parameter which affect the estimation accuracy and speed, because within one data sample, the power density is calculated simultaneously, which has been introduced and tested in Section 2.4.3. 1 u(i)uT (i − ) l

(10)

1 y(k)uT (i − ) l

(11)

l−1

Ruu () =

i=1 ∞

Ryu () =

of the unexcited control signals will be used to validate the system identification model. 3.2. Building energy system identification model performance index Model forecasting accuracy when compared with EnergyPlus results and speed are two most important indices in this study. Coefficient of determination, R2 (Eq. (12)) is used to measure the forecasting data accuracy [36]. R2 =

=0

Fig. 7 shows the procedure of using system identification model (SID model hereafter) to develop the on-line building energy model. During the training process, the excited building control signals will be generated in Matlab according to the exciting scheme discussed above and sent to EnergyPlus through BCVTB. The EnergyPlus simulation results of the excited system will be used to train the SID model, and calculate its Markov parameters for each input in transfer function. On the other hand, the EnergyPlus simulation results

n ¯ yi − yˆ¯ ) (yi − y)(ˆ n i=1 2 n ¯ (y − y) (ˆy − yˆ¯ ) i=1 i i=1 i

(12)

where yi is the energy consumption data from EnergyPlus, yˆ i is the energy consumption forecasting data, y¯ and yˆ¯ are their average. 4. System identification model estimation results Two different models, core zone model and perimeter zone model, have been developed and applied to model and forecast the

Fig. 6. System identification model development procedure.

X. Li, J. Wen / Energy and Buildings 82 (2014) 1–12

7

Fig. 7. On-line building energy model development procedure.

heating/cooling load and energy, individually. The training period for building load and heating energy estimation models is from January 1st to January 10th, and the forecasting period for them is from January 11th to January 13th. While the training period for cooling energy estimation models is July 1st to 10th, and the forecasting period is July 11th to 13th. During training period, all the system exciting signals are applied into the building model, but regular control signals are used in the building during the forecasting/testing period (Fig. 8). 4.1. Building heating and cooling load off-line estimation results In order to test the SID model’s capability in capture the thermal response dynamics of building thermal mass, especially the solar heat gains, an off-line building heating and cooling load estimation model is developed. This testing study is based on the EnergyPlus simulation for from January 1st to January 10th. As discussed in pervious sections, for both core zone and perimeter zone models the excitation signals are injected into the building model every 30 min. The time step of training data generation EnergyPlus model is 5 min, while that of SID model is 15 min. The shorter training model time step is to better capture the building dynamics. With considering the dynamics difference between core zone and perimeter zone, the data sampling length are 240 and 180 min, respectively. Since building heating/cooling load is hard to measure in real field, the building load estimation is an off-line process. The load estimation

building is trained based on the EnergyPlus simulation results for January 1st to January 10th. In Fig. 9a and b, ten days’ heating and cooling load training results from the system identification model are compared with those from EnergyPlus model. The results in this figure are the amount of heating load subtracting that of cooling load. Hence the value above zero is heating load, while the value below zero is cooling load. Since the internal heat gains are very large in the building, the core zone just need some heating in the morning. After all the equipment turning on, core zone requires cooling in the afternoon. Similar phenomenon is identified for the perimeter zones, but in the morning more heating is requested. From this plot, we can see the green line (SID) mostly capture dynamics of the blue line (EnergyPlus). The seventh day of the training is a holiday, when all the HVAC system is off. For the simplification of the presentation, only the results of core zone and west perimeter zone are illustrated here. Those estimation results of other perimeter zones are very similar to that of west perimeter zone. Without considering the outliers, the R2 value of the training data estimation is summarized in Table 4. Both of these two are above 0.95. Once the building load system identification model is trained to have good estimation accuracy, it is then used to forecast the future building load. Fig. 9c plots the building load forecasting results using the system identification model discussed above. In the perimeter zone load estimation result plot of Fig. 9d the EnergyPlus simulated load (Qep) is zero at night when system if off. However, forecasted load from SID model (Qes) does not have a system on/off identifying variable. It has been added into cooling energy estimation model. Therefore, SID model load estimation model cannot identify the system off period. It will calculate the building load as the system is on, which the energy estimation model will read in the system on and off schedule first and then calculate the energy consumption. In the perimeter zone simulation results plots, we can see that our model will overestimate the cooling load and energy at noon from 11 am to 2 pm and overestimate the heating load and energy at late afternoon from 4 pm to 10 pm. That is because of this on-line model does not consider the

Table 4 Building load system identification model training estimation accuracy.

Fig. 8. Temperature setpoints and equipment schedule during forecasting period.

R2

Training

Forecasting

Core zone Perimeter zone

0.937 0.961

0.953 0.948

8

X. Li, J. Wen / Energy and Buildings 82 (2014) 1–12

Fig. 9. Building heating and cooling load estimation results: (a and b) Building heating/cooling load training results; (c and d) Building heating/cooling load fore casting results; (e and f) Building load forecasting error.

heat storage in the building thermal mass, which will cause the lag effect for building heating and cooling load. And Fig. 9e and f plot the discrepancy between the results from the system identification model (Qes) and those from EnergyPlus simulation model (Qep), all the blue circles are clustered around the line: Qes = Qep. From these two figures, we can see the forecasting accuracy of core zone is better than that of perimeter zone, which coordinate to the load estimation. And the on-line model will overestimate the load at 8 am when the system is starting up. Those are the outliers in Fig. 9e and f, which are acceptable, because the system has a very huge dynamic change at the starting up and shutting down periods. And this period is very short, which will not affect the overall results too much.

4.2. Building cooling energy on-line estimation results Different to heating/cooling load, building cooling energy consumption is easy to measure in real field. The heating/cooling energy forecasting model is essential to building model predictive control for energy saving. Heating system coefficient is a constant value, so heating energy consumption is very easy to estimate based on the heating load estimation. It is not necessary to develop a new SID model for heating energy forecasting. Therefore, this study will only focus on the cooling energy modeling and forecasting. One very important variable was added into energy estimation model, which is supply fan heat into the air stream. Fig. 10 shows the plots of cooling energy forecasting results. Similar to the results

X. Li, J. Wen / Energy and Buildings 82 (2014) 1–12

9

Fig. 10. Building cooling energy forecasting results: (a) East zone cooling energy forecasting; (b) South zone cooling energy forecasting; (c) West zone cooling energy forecasting; (d) North zone cooling energy forecasting; (e) Core zone cooling energy forecasting.

of load estimation, the model would also overestimate the energy consumption when the system is starting up and shutting down. The accuracy of heating energy forecasting results is better than that of cooling energy forecasting (Table 5). That is because the heating energy efficiency coefficient is usually a constant value, while the cooling energy efficiency is depended on the partial load ratio, and zone temperature setpoint, etc. The core zone cooling energy forecasting results are better than training results, because during training period, the system has been excited, which has much more dynamics than that at forecasting period. These five plots show that the on-line model

Table 5 System identification model forecasting accuracy. R2

East zone South zone West zone North zone Core zone Whole building

Accuracy (R2 )

Speed (S)

Training

Forecasting

Training

Forecasting

0.967 0.971 0.969 0.954 0.965 0.944

0.957 0.963 0.966 0.943 0.955 0.956

18.30 19.54 19.90 18.93 18.64 47.33

0.0041 0.0094 0.0042 0.0043 0.0035 0.0089

10

X. Li, J. Wen / Energy and Buildings 82 (2014) 1–12

Fig. 11. Whole building cooling energy forecasting results.

underestimates the cooling energy consumption, especially in west zone and zone, at afternoon when direct solar radiation is strong. This is because solar radiation at afternoon is very strong and the solar angle is keep changing which brings in huge dynamics of the building energy system. However, east perimeter zone model overestimates the cooling energy because the solar only directly radiate on the east envelop at morning, and direct solar radiation at morning is very small. In additional to the direct solar radiation changing, the building thermal mass heat storage also affects the cooling energy consumption. The solar radiation transmitted into the building would not become cooling load of the HVAC system immediately. It will be first absorbed by the internal walls, floor, furniture and etc., and then released out by convection and radiation. 4.3. Whole building energy estimation results The models discussed above are zone based models, which forecast building cooling energy consumption zone by zone. Since the building studied in this paper is a one-floor small building, an overall model is also developed for the whole building. In this whole building model, the system identification method is the same as that of zone model, including system excitation signals, frequency response function model, just with model input and output adjust. The output is the cooling energy of the whole building, and the inputs are all the inputs of core zone model and perimeter zone models. The whole building energy consumption forecasting results from the whole building model are plotted in Fig. 11a and b. Fig. 11b illustrates the comparison of EnergyPlus simulated whole building energy consumption (Qes) and SID model forecasted while building energy consumption (Qep). Due to the underestimation of direct solar radiation related cooling energy consumption at afternoon when cooling load is very high, a large discrepancy between Qes and Qep exists, as highlighted in green dashed circles. Even with this underestimation, the overall forecasting accuracy is still desirable (Table 5). Simulation speed is another crucial factor for MPC in practical use. The training and forecasting speed of SID model are summarized in Table 5. The training calculation time is relative longer because the power density is evaluated step by step. The forecasting calculation time is trivial, because the forecasting model is linear equation calculation of input data and Markov parameters saved from training period. Fortunately, if this model is applied in a MPC model, the training process just need to be conducted once and all the Markov parameters and transfer functions would be saved

and stay unchanging. Therefore, the forecasting period model will run multiple times as needed in MPC model. Therefore, this model is still suitable for any “searching” based optimization and control algorithms. 5. System identification model application in medium office building In order to test the accuracy of the system identification modeling strategy, the approach developed and validated in previous sections has been utilized in another medium size commercial building. The building structure and mechanical system are different to the small building tested above, which is introduced in detail in following sections. 5.1. Medium office building description The system identification modeling approach developed in this study has been utilized in another medium office building, in order to test its accuracy and robustness. This medium building is a 4982.2 m2 (53,628 ft2 ), three-story, 15-zone building [24]. The window-to-wall ratio of this building’s facades is approximately 33.0%, and the windows are also equally distributed. The thermal properties of the windows of this building are the same as those in the previous small commercial building. The mechanical configurations of this building are tabulated in Table 6. Different to the small commercial building, this building is using variable-air-volume (VAV) AHU. Primary cooling is also provided by electricity through DX coils. The coefficients of performance (COP) of these cooling coils are 3. Primary heating is provided by a natural gas (NG) boiler, which have a70% annual fuel efficiency utilization (AFUE). Zone reheat is provided electric resistance heating.

Table 6 Mechanical systems of the medium building. Specifications System Main cool coil Main heat coil Zone reheat Heat plant Heat efficiency

3 VAV, AHUs DX, COP 3 NG Furnace Electric Gas Central Boiler 70% AFUE

X. Li, J. Wen / Energy and Buildings 82 (2014) 1–12

11

Fig. 12. Medium building cooling energy estimation: (a) First floor model training results; (b) First floor cooling energy forecasting; (c) Second floor cooling energy forecasting; (d) Whole building cooling energy forecasting.

5.2. Cooling energy on-line forecasting using system identification model The system identification model developed and validated in previous sections have been applied in this building with input adjust based on the system configurations. Since there are 15 thermal zones in this building, it is very time consuming to develop 15 different models to forecast the cooling energy consumption. Therefore three different models were developed for each floor which consists of five zones: core zone, east perimeter zone, south perimeter zone, west perimeter zone and north perimeter zone. Same system excitation strategy was used in this building for building temperature setpoints and indoor equipment schedules. The inputs and outputs of each model are also similar to the small building model. But in each model, temperature of each room and each adjacent room at each floor are included.

5.3. Cooling energy forecasting results As introduced in Section 5.2, building cooling energy consumption was forecasted at each floor individually. Fig. 12 shows the training and forecasting results. The forecasting accuracy and speed are summarized in Table 7, which illustrates that the forecasting results capture the trend of real energy consumption measurement. The accuracy and speed are still acceptable for MPC, even

though they are not as good as those for the small building, because medium building has more disturbances and it requires longer calculation time to include all these disturbances into the model. It is illustrated in Fig. 12 that the on-line estimation model underestimates the cooling energy consumption during the HVAC starting up and shutting down periods. Comparing to the forecasting results of the small building, the forecasting results in this medium building has lagers errors. That is because there are more disturbances in the medium building than in small buildings. However, the on-line model still capture the overall trend of the real energy consumption and the accuracy is still above 88%. Next step of this research will focus on updating forecasting results based on real measurements when it is necessary to improve the accuracy, especially at the high dynamics period.

Table 7 System identification model forecasting accuracy and speed for medium building. Accuracy (R2 )

First floor Second floor Third floor Whole building

Speed (S)

Training

Forecasting

Training

Forecasting

0.931 0.941 0.951 0.949

0.878 0.871 0.877 0.872

42.6 52.9 42.9 97.3

0.011 0.015 0.013 0.026

12

X. Li, J. Wen / Energy and Buildings 82 (2014) 1–12

6. Conclusion and future work This study introduced a novel systematic methodology for online building energy estimation model development and validation. A system excitement scheme and a Matlab-BCVTB-EnergyPlus testbed was developed and validated for system identification model development. The excitement scheme is able to guarantee enough data at the high frequency and low frequency around the building operation range, which can be applied in any other system identification models for building energy simulation. Frequency response function model with realized by power spectral density model was implanted to forecasting building load and energy consumption. This on-line building energy model can achieve over 95% forecasting accuracy within 1 s in a small building case, which is suitable for any on-line building operation and MPC model. This modeling approach has also been applied in a medium office building. Limited to the calculation speed, the cooling energy forecasting accuracy is still around 88%. We expect to further improve the accuracy of this on-line estimation model by introducing a thermal mass storage variable to simulate the thermal storage effect and a solar angle parameter to estimate effect of direct solar radiation on energy consumption. Data fusion techniques, such as Kalman filter, is planned to further improve the forecasting accuracy, reduce the calculation time and improve the model robustness. In this study, only DX coil with CAV AHUs have been modeled, but we plan to model other mechanical system such as VAV and gas heating. In additional, this on-line building model is planned to be applied in some optimization models to determine the optimal building operation strategies. Acknowledgement Financial support provided by the U.S. National Science Foundation Award 1239247 is greatly appreciated. References [1] U.S. DOE, Buildings Energy Data Book, http://buildingsdatabook.eren.doe.gov/, 2011. [2] S. Katipamula, M.R. Brambley, Review article: methods for fault detection, diagnostics, and prognostics for building systems—a review, Part I, HVAC&R Research 11 (1) (2005) 3–25. [3] National Energy Technology Laboratory (NETL), Demand Dispatch–Intelligent Demand for a More Efficient Grid, in: DOE/NETL-DE-FE0004001, Morgantown, West Virginia, August 2011. [4] U.S.D.o.E. Office of Energy Efficiency and Renewable Energy, EnergyPlus Energy Simulation Software: About EnergyPlus. [5] M. Wetter, P. Haves, A Modular building controls virtual test Bed for the Integration of heterogeneous systems, in: Proceedings of 3rd Nation Conference of IBPSA-USA SimBuild, Berkeley, California, 2008. [6] J. Ma, J. Qin, T. Salsbury, P. Xu, Demand reduction in building energy systems based on economic model predictive control, Chemical Engineering Science 67 (1) (2012) 92–100. [7] C.D. Corbin, G.P. Henze, P. May-Ostendorp, A model predictive control optimization environment for real-time commercial building application, Journal of Building Performance Simulation 6 (3) (2013). [8] E.E.a.R.E.O., U.S. Department of Energy, Building Technology Program, Net-Zero Energy Commercial Building Initiative, Commercial Building Benchmark Models, 2009, Available from: http://www1.eere.energy.gov/buildings/ commercial initiative/benchmark models.html [9] M. Wetter, Design optimization with GenOptR , in: Building Energy Simulation, Lawrence Berkeley National Laboratory, 2000, September/October.

[10] B. Coffey, F. Haghighat, E. Morofsky, E. Kutrowski, A software framework for model predictive control with GenOpt, Energy and Buildings 42 (7) (2010) 1084–1092. [11] W.J. Cole, E.T. Hale, T.F. Edgar, Building energy model reduction for model predictive control using OpenStudio, in: 2013 American Control Conference, Washington, DC, 2013. [12] K. Yun, R. Luck, P.J. Mago, H. Cho, Building hourly thermal load prediction using an indexed ARX model, Energy and Buildings 54 (2012) 225–233. [13] J. Yang, H. Rivard, R. Zmeureanu, On-line building energy prediction using adaptive artificial neural networks, Energy and Buildings 37 (12) (2005) 1250–1259. [14] Y. Chen, P.B. Luh, C. Guan, Y. Zhao, L.D. Michel, M.A. Coolbeth, P.B. Friedland, S.J. Rourke, Short-term load forecasting: similar day-based wavelet neural networks, IEEE Transactions on Power Systems 25 (1) (2010) 322–330. [15] J.E. Braun, N. Chaturvedi, An inverse gray-box model for transient building load prediction, HVAC&R Research 8 (1) (2002) 73–99. [16] F. Oldewurtel, A. Parisio, C.N. Jones, D. Gyalistras, M. Gwerder, V. Stauch, B. Lehmann, M. Morari, Use of model predictive control and weather forecasts for energy efficient building climate control, Energy and Buildings 45 (2012) 15–27. [17] K.-H. Lee, J.E. Braun, Development of methods for determining demandlimiting setpoint trajectories in buildings using short-term measurements, Building and Environment 43 (10) (2008) 1755–1768. [18] S. Wang, X. Xu, Simplified building model for transient thermal performance estimation using GA-based parameter identification, International Journal of Thermal Sciences 45 (4) (2006) 419–432. [19] D. Kim, J.E. Braun, Model-based predictive control for buildings with decoupling and reduced-order modeling, in: 2nd International High Performance Buildings Conference, Purdue University, 2012. [20] J. Ma, S.J. Qin, T. Salsbury, Model predictive control of building energy systems with balanced model reduction, American Control Conference (ACC), 2012, IEEE (2012) 3681–3686. [21] Z. O’Neill, S. Narayanan, R. Brahme, Model-based thermal load estimation in buildings, in: Proceedings of SimBuild, 2010. [22] L. Ljung, System Identification, Wiley Online Library, 1999. [23] S. Privara, Z. Vána, D. Gyalistras, J. Cigler, C. Sagerschnig, M. Morari, L. Ferkl, Modeling and identification of a large multi-zone office building, IEEE International Conference on Control Applications (CCA), IEEE (2011) 55–60. [24] C. Agbi, Z. Song, B. Krogh, Parameter identifiability for multi-zone building models, IEEE 51st Annual Conference on Decision and Control (CDC), IEEE (2012) 6951–6956. [25] M. Deru, K. Field, D. Studer, K. Benne, B. Griffith, P. Torcellini, B. Liu, M. Halverson, D. Winiarski, M. Rosenberg, US Department of Energy Commercial Reference Building Models of the National Building Stock, 2011. [26] S.P. Corgnati, E. Fabrizio, M. Filippi, V. Monetti, Reference buildings for cost optimal analysis: method of definition and application, Applied Energy 102 (2013) 983–993. [27] L. Hendricken, K. Otto, J. Wen, P. Gurian, W. Sisson, Capital costs and energy savings achieved by energy conservation measures for office buildings in the greater Philadephia Region, in: SimBuild Madison, WI, 2012. [28] L.C. Ng, A. Musser, A.K. Persily, S.J. Emmerich, Indoor air quality analyses of commercial reference buildings, Building and Environment 58 (2012) 179–187. [29] F. Zhao, I.J. Martinez-Moyano, G. Augenbroe, Agent-based modeling of commercial building stocks for policy support, in: The 12th International Building Performance Simulation Association Conference, Sydney, Australia, 2011. ˇ E. Zˇ áˇceková, J. Cigler, Building modeling: selection of the [30] S. Prívara, Z. Vána, most appropriate model for predictive control, Energy and Buildings 55 (2012) 341–350. [31] J.N. Juang, Applied System Identification, Prentice Hall, Englewood Cliffs, New Jersey, 1994. [32] M. Braun, R. Ortiz-Mojica, D. Rivera, Application of minimum crest factor multisinusoidal signals for “plant-friendly” identification of nonlinear process systems, Control Engineering Practice 10 (3) (2002) 301–313. [33] M. Braun, D. Rivera, A. Stenman, W. Foslien, C. Hrenya, Multi-level pseudorandom signal design and “model-on-demand” estimation applied to nonlinear identification of a RTP wafer reactor, in: Proceedings of the American Control Conference 1999, IEEE, 1999, pp. 1573–1577. [34] D.E. Rivera, H. Lee, H.D. Mittelmann, M.W. Braun, Constrained multisine input signals for plant-friendly identification of chemical process systems, Journal of Process Control 19 (4) (2009) 623–635. [35] D.E. Rivera, X. Chen, D.S. Bayard, Experimental design for robust process control using schroeder-phased input signals, American Control Conference, 1993, IEEE (1993) 895–899. [36] J.L. Devore, Probability and Statistics for Engineering and the Sciences, Brooks/Cole Publishing Co., 2004.