The role of buildings in a changing environmental ERA: a modeling approach to quantifying and reducing missions

As the scientific community continues to make predictions about the implications of a changing climate within our social and material infrastructures, the call for solutions that mitigate these implications is becoming more urgent than ever before. Responsible for a large portion of the energy demand in developed regions, buildings are central to the problem and the solution. In this work, a tool is developed for quantifying indirect emissions and water usage associated with power drawn from the electrical grid. The tool uses the WattTime API and databases published by the U.S. Environmental Protection Agency and the U.S. Department of Energy. It is developed as an open-access tool that provides hourly externality data for select locations within the United States. This kind of data can be hugely influential in understating the impact of building operations and the extent to which this extends beyond the immediate surroundings of a building. After a brief focus on the environmental issues that buildings can pose, the thesis moves into an example of how Model Predictive Control (MPC), a control methodology gaining popularity in building energy control, can be used to mitigate the burden that buildings can have on grid operations and the environment. Building energy modeling is the focal point of MPC algorithms, and for many reasons, is the deciding factor in the feasibility of MPC in real building systems. In particular, uncertainty in complex building models and in the data that are used to make them act as a barrier for wide-scale infiltration into the building sector. Here, an uncertainty analysis for a small office building is used to better understand the extent to which models can be simplified without compromising accuracy and how uncertainty in the most important building modeling parameters can impact the MPC optimization procedure. It is intended that this work is presented through the lens of an applied framework and that the tools and methods herein extend our opportunities for understanding and mitigating the environmental challenges that we are experiencing.

THE ROLE OF BUILDINGS IN A CHANGING ENVIRONMENTAL ERA: A MODELING APPROACH TO QUANTIFYING AND REDUCING EMISSIONS by Kaden Ezra Plewe A thesis submitted to the faculty of The University of Utah in partial fulfillment of the requirements for the degree of Master of Science Department of Mechanical Engineering The University of Utah August 2019 Copyright c Kaden Ezra Plewe 2019 All Rights Reserved The University of Utah Graduate School STATEMENT OF THESIS APPROVAL The thesis of Kaden Ezra Plewe has been approved by the following supervisory committee members: Amanda D. Smith , Chair(s) 30 May 2019 Date Approved Kent S. Udell , Member 17 May 2019 Date Approved Mingxi Liu , Member 17 May 2019 Date Approved by Bruce K. Gale , Chair/Dean of the Department/College/School of Mechanical Engineering and by David B. Kieda , Dean of The Graduate School. ABSTRACT As the scientific community continues to make predictions about the implications of a changing climate within our social and material infrastructures, the call for solutions that mitigate these implications is becoming more urgent than ever before. Responsible for a large portion of the energy demand in developed regions, buildings are central to the problem and the solution. In this work, a tool is developed for quantifying indirect emissions and water usage associated with power drawn from the electrical grid. The tool uses the WattTime API and databases published by the U.S. Environmental Protection Agency and the U.S. Department of Energy. It is developed as an open-access tool that provides hourly externality data for select locations within the United States. This kind of data can be hugely influential in understating the impact of building operations and the extent to which this extends beyond the immediate surroundings of a building. After a brief focus on the environmental issues that buildings can pose, the thesis moves into an example of how Model Predictive Control (MPC), a control methodology gaining popularity in building energy control, can be used to mitigate the burden that buildings can have on grid operations and the environment. Building energy modeling is the focal point of MPC algorithms, and for many reasons, is the deciding factor in the feasibility of MPC in real building systems. In particular, uncertainty in complex building models and in the data that are used to make them act as a barrier for wide-scale infiltration into the building sector. Here, an uncertainty analysis for a small office building is used to better understand the extent to which models can be simplified without compromising accuracy and how uncertainty in the most important building modeling parameters can impact the MPC optimization procedure. It is intended that this work is presented through the lens of an applied framework and that the tools and methods herein extend our opportunities for understanding and mitigating the environmental challenges that we are experiencing. For the generation that doesn’t have a voice. CONTENTS ABSTRACT . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii LIST OF TABLES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii LIST OF FIGURES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ix ACKNOWLEDGEMENTS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi 1. 2. 3. INTRODUCTION . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 1.1 Externalities in the Electric Power Industry . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.1 Emissions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.1.2 Water Usage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2 Buildings . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.1 Building Performance Simulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.2.2 Building Controls . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1.3 Materials and Resources . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1 2 2 3 3 4 5 QUANTIFYING EXTERNALITIES IN THE POWER SECTOR . . . . . . . . . . . . . . 6 2.1 Electrical Grid Operation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 Power Externality Correlation Tool . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.1 Materials and Methods . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.2 Software Architecture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2.3 Software Functionalities . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Case Studies . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.4 Implications . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7 11 11 14 16 16 19 20 MODEL PREDICTIVE CONTROL IN BUILDING THERMAL SYSTEMS . . . . 21 3.1 Problem Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 System Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2.1 Thermal RC Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.1 Basic On/Off Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3.2 Model Predictive Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.1 Comparison to Basic On/Off Control . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.4.2 Demand Peak Shaving . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4. 22 23 24 28 30 31 35 36 41 43 UNCERTAINTY IN BUILDING PERFORMANCE SIMULATION . . . . . . . . . . . 45 4.1 Problem Definition . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1 System Model . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1.1 Nominal Performance . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.1.1.2 Input Parameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Uncertainty and Sensitivity Calculation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Building Performance Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4 Results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.1 Uncertainty Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.2 Sensitivity Analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.4.3 Uncertainty and Sensitivity Analysis Comparison . . . . . . . . . . . . . . . . . . 4.4.4 Uncertainty Impact on Optimization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5. 47 48 49 50 51 54 56 57 59 62 65 68 CONCLUDING REMARKS . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70 REFERENCES . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72 vi LIST OF TABLES 2.1 Generation Data Availability from WattTime API [1] . . . . . . . . . . . . . . . . . . . . . . 12 2.2 Water Consumption and Withdrawal Factors [2] . . . . . . . . . . . . . . . . . . . . . . . . . 13 2.3 Calculated emission Factors for CO2 , NOx and SO2 using Equation 2.3. Emission factors are calculated with eGRID data [3] for Coal, Gas and Biomass for each balancing authority included in the WattTime API [1] . . . . . . . . . . . . . . . . 14 3.1 Heating design power calculated by EnergyPlus for each thermally regulated zone for a design temperature of -16.30 ◦ C. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2 Data used to calculate the thermal capacitance of each zone. . . . . . . . . . . . . . . . 28 3.3 Total thermal resistance values obtained using resistance and areas listed in the EnergyPlus small office model. Values are oriented to represent the resistance between the row and column nodes. . . . . . . . . . . . . . . . . . . . . . . . . . . 28 3.4 Heating and cooling setpoint schedules used over the winter control simulation period. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 29 3.5 Constant heat supply levels that are used for on/off control schemes. . . . . . . . . 30 3.6 Basic on/off control algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 31 3.7 Model Predictive Control algorithm used for on/off and variable heating supply MPC simulations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 33 3.8 Total heating energy usage for 48-hour period in January five different control simulations. MPC10, MPC20 and MPC30 are MPC variable control simulations using a forecasting horizon of 10, 20 and 30, respectively. . . . . . . . . . . . . . 36 4.1 Building performance simulation advancement timeline for the 1960s through present day. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46 4.2 Predicted mean vote sensation scale used for thermal comfort quantification. . 49 4.3 Nominal building performance metrics calculated over two-weeks in the summer and winter. Two-week sum (facility electric energy) and two-week average (predicted mean vote). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50 4.4 EnergyPlus IDF classes included in UA/SA. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.5 Brief of algorithm used to apply simulated uncertainty intervals to input parameters. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52 4.6 Uncertainty analysis for primary energy, emissions and cost quantification over a two-week summer simulation period. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 4.7 Statistical comparison of output distributions for EnergyPlus model and Gaussian process model for 4000 training samples. . . . . . . . . . . . . . . . . . . . . . . . 61 4.8 Top 10 parameters, based on DGSM index, for facility electric energy output. . 62 4.9 Top 10 parameters, based on DGSM index, for predicted mean vote output. . . 63 4.10 Statistical comparison between EnergyPlus simulations with 20% uncertainty applied to all parameters and EnergyPlus simulations with 20% uncertainty applied to the top 10 parameters and top 5 parameters for FEE and PMV using a sample size of 5000. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 4.11 Statistical comparison of output distributions for EnergyPlus model and Gaussian Process Regression (GPR) model for 4000 training samples using only the top 10 and top 5 parameters identified using DGSM sensitivity indices (4.8, 4.9). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 64 viii LIST OF FIGURES 2.1 Map of NERC interconnections and regional entities [4] . . . . . . . . . . . . . . . . . . . 8 2.2 Map of the subregions represented in the latest eGRID database published in 2014. Four subregions are excluded in this image: AKMS and AKGD, in Alaska and HIMS and HIOA, in the Hawaiian Islands. Crosshatching indicates that the geographical area covers multiple eGRID regions [3]. . . . . . . 9 2.3 PECT architecture outlining the process that is used to obtain emissions and water usage estimates. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15 2.4 Total generation, emissions and water usage attributed to various fuel types reported by the MISO balancing authority. Software query was performed for Kalamazoo, Michigan for April 22, 2018. The reg. fit lines are obtained using relaxed least squares regression. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17 2.5 Generation mix for three locations: (a) Michigan (top), (b) Virginia (bottom left), and (c) Vermont (bottom right) on July 17, 2017. The balancing authorities for these locations are (a) MISO, (b) PJM, and (c) ISONE. . . . . . . . . . . . . . . . . . . 18 3.1 Isometric and top view of the interior of the small office prototype building used in this study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2 Thermal RC circuits for external and core zones with heating input and the attic zone. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25 3.3 Outdoor dry bulb temperature for a 48-hour period in January with the thermal response of the zone temperatures without heat input. . . . . . . . . . . . . . 29 3.4 Control input and state response trajectory for basic on/off control over a 48-hour period (a) and for a one-hour period (b) in January. . . . . . . . . . . . . . . . . 38 3.5 Control input and state response trajectory for MPC on/off control over a 48-hour period (a) and for a one-hour period (b) in January. . . . . . . . . . . . . . . . . 39 3.6 Control input and state response trajectory for MPC variable control over a 48-hour period (a) and for a one-hour period (b) in January. . . . . . . . . . . . . . . . . 40 3.7 Forecast horizon comparison for MPC variable control simulation with a forecast horizon of 10, 20 and 30 minutes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41 3.8 Comparison of temperature trajectories for MPC control utilizing a difference cost weighted with φ = 0, φ = 1000 and φ = 10000. . . . . . . . . . . . . . . . . . . . . . . 42 3.9 Comparison of heat supply trajectories for MPC control utilizing a difference cost weighted with φ = 0, φ = 1000 and φ = 10000. . . . . . . . . . . . . . . . . . . . . . . 43 4.1 Isometric (left) and top view of the interior (right) of the small office prototype building used in this study. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 48 4.2 Nominal performance for small office building simulations over summer (07/16-07/30) and winter (01/15-01/29) simulation periods. . . . . . . . . . . . . . . . 50 4.3 Interpretation of box plots that are used to show output distributions. . . . . . . . 56 4.4 Daily trends and 95% confidence interval relating to each category that 20% simulated uncertainty was applied to separately. . . . . . . . . . . . . . . . . . . . . . . . . . 57 4.5 Output distributions for simulated 20% uncertainty the three separate categories and for results related to simulating uncertainty in all parameters. . . . . . 58 4.6 Gaussian process model convergence analysis for two-week summer simulation. Bhattacharyya distance as a function of training sample size is plotted to represent the convergence of EnergyPlus (EP) and Gaussian Process (GP) model output distributions for 5000 samples. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.7 Comparison of EnergyPlus model with Gaussian process model that was trained with 4000 simulation samples. Output distributions are for a sample size of 5000 points. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60 4.8 Derivative-based Global Sensitivity indices for facility electric energy output and predicted mean vote output. All indices are included in the left two subplots and the top 10 parameters for both output variables are included in the right two subplots. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 61 4.9 Comparison of EnergyPlus model with Gaussian process models that were trained with all 907 parameters, top 10 parameters and top 5 parameters based on the DGSM index ranking. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63 4.10 Comparison of output uncertainty for EnergyPlus simulations with simulated input uncertainty in three input categories, top 10 parameters and top 5 parameters based on the DGSM indices. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65 4.11 Comparison of Pareto frontiers achieved with the optimization procedure repeated four times for five different uncertainty scenarios: no uncertainty and 20% uncertainty applied to the top 5 and 10 parameters for the PMV and the FEE output. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67 x ACKNOWLEDGEMENTS I would like to recognize the people in my life who have supported me throughout my time at the University of Utah. My work here and the many experiences that are behind it would not have been possible without them. Firstly, I want to express appreciation for my advisor, Dr. Amanda D. Smith. She has played a critical role in giving me the tools, guidance and confidence necessary for carrying out this research and I hope that I can utilize the things I have learned from her to have a substantial impact in my future work. I would also like to thank my committee members Dr. Mingxi Liu and Dr. Kent S. Udell, as their instruction and assistance near the end of the completion of this work was critical. The Global Change and Sustainability Center and Northrop Grumman together have given me the financial freedom to pursue topics that I am passionate about, and I thank them for that. Thank you to my siblings, Raegan, Reagan, Elliott, Cohen and Isabella for your support and understanding with regard to the lost time spent together as a result of my studies. Thank you Mom, Erica, Dad and Jillian for providing me with the love, support and opportunities that have brought me to this point in life. And lastly, thank you to Jermy Thomas, Zahra Fallahi, Aowabin Rhaman, Zahra Ghaemi, Jaron Peck, Carlo Bianchi and all of the friends and colleagues who aren’t mentioned here explicitly. You have been a tremendous help at every step and have made it all enjoyable. CHAPTER 1 INTRODUCTION This thesis seeks to pose the building energy sector as an artifact of modern infrastructure that is an interconnected component of global problems and solutions. With the changes in climate and resource availability that we expect to see in the future, building energy systems will need to be designed and retrofitted in a way that meets comfort expectations while working within certain limitations. This is the interaction between buildings and the environment that is being highlighted in the title of this thesis. With evolving building energy codes, the adoption of intermittent renewable sources, and the potential adoption of carbon pricing, building energy systems will need to become more efficient, interface with demand management and real-time demand response programs and be generally improved in terms of data availability and control flexibility. The collective advancement of this field, to accomplish environmental, economic, and social objectives, will require an interdisciplinary approach. This chapter provides a brief overview of the interconnected topics that are discussed throughout the body of the thesis. 1.1 Externalities in the Electric Power Industry The electric power industry is the backbone of our modernized world and has the task of providing reliable electricity to a modular demand base in a time of use framework. In the United States, the supply of electricity to individual users is accomplished through a connection between end users and the main supply system that involves numerous points of management. A more in-depth description of how this is managed is covered in Chapter 2. Power providers that supply electricity to a regional grid make up the regional generation portfolio and can be found in publicly accessible databases [3]. At any point in time, only a fraction of the generation portfolio is actively supplying electricity to the grid, and the makeup of different types of power plants that are active at any moment is defined as the fuel mix. This introduces a temporal dependency to the cost and externalities of grid 2 electricity because every power plant has unique attributes that affect these things such as the fuel they use, the type of energy conversion cycles they use, their age, ramp up cost and location. There is also a spatial variability associated with electricity purchased from the grid that is based on the site’s physical distance to active power providers. This is a result of both the physical constraints associated with the transfer of electrons and directing capabilities of grid managers. It is important to recognize this temporal and spatial variation because discussions involving cost and impacts associated with the electric power industry can vary greatly based on these temporal and spatial variabilities. 1.1.1 Emissions The emissions associated with the electric power industry are commonly discussed in terms of carbon dioxide. This is because CO2 makes up the largest fraction of the total emissions associated with combustion processes and is a well-known greenhouse gas. Nitrous oxides and Sulfur dioxide emissions are also typically included in power plant emission reporting. The electric power industry is responsible for approximately 38% of CO2 emissions, 13% of NOx emissions and 59% of SO2 emissions in the United States [5]. The emissions associated with power plants are considered indirect for any electric end user, and they are considered direct for the power plant itself. Because they are hidden at the point of use, indirect emissions are often disregarded by end users. It is essential that indirect emissions are considered for buildings because they make up the largest portion of total building life cycle emissions [6]. 1.1.2 Water Usage Water is often an overlooked commodity in the evaluation of the electric power industry. But it has a key role in determining the operational feasibility of a power plant in a given region because the flow rates that are required in the cooling process are far larger than what can be provided by utility water pipes. For this reason, power plants often take advantage of natural water resources such as rivers, lakes, or oceans for their cooling needs. The United States Geological Survey estimates that approximately 45% of total withdrawals are associated with the power sector [7]. Water withdrawals refer to the total amount of water that is withdrawn from its source. Water consumption is the portion of water withdrawn that isn’t returned to its source, for which the power sector accounts for a much smaller 3 percentage (<10%). These withdrawals mainly come from thermoelectric power plants that use water in the condensing process of the power cycle, which either operate by passing water through the condenser and then returning it directly to the source or by recycling water repeatedly. The water sector and the power sector have a unique relationship that is sometimes referred to as the water-energy nexus. The water-energy nexus describes an interdependency between the two industries, where the water supply system has a high energy demand, and a large portion of the power generation portfolio has a high water demand. This interdependency is highly concerning when we consider the implications of expected water scarcity, especially in places like the Southwest United States [8]. 1.2 Buildings As an essential component of modern human activity, buildings are perhaps one of the most interconnected aspects of our culture. The context that buildings provide in our daily lives has implications in social class (size, amenities), popular culture (architecture, style), religion (temples, mosques) and the comfort and productivity that we seek to maintain. But beyond this, buildings have a high energy demand and place a heavy burden on our environment. In developed countries, as much as 40% of total final energy consumption is associated with buildings [9], [10]. In the past century, the community of building researchers and practitioners has made significant advancements in building energy efficiency, thermal comfort, building design, and building controls. As a result, we benefit from standards that ensure our comfort and health when we occupy closed spaces. Advancements in building science include high-power computing, machine learning, and automation to optimize the performance of buildings with respect to many objectives (thermal comfort, cost, emissions). These advancements are also being applied to take advantage of interconnected microgrids and distribution networks, where demand-side management is utilized to balance loads and reduce total energy usage. 1.2.1 Building Performance Simulation Building performance simulation (BPS), in the context of building energy systems, seeks to capture the behavior of a building so that its energy performance can be predicted based on expected temperature and occupancy disturbances. BPS has traditionally been 4 used in the initial design and commissioning phase of a building life cycle. But in the past decade, it has gained relevance in the continuous operation and control of buildings due to the development of model-based control techniques for building energy systems. The complexity of a building model can vary based on the application and requirements associated with the end use of the building model. There is a large selection of software tools that are designed to assist the process of BPS [11], with EnergyPlus [12] being one of the most commonly utilized. To accurately model a building, one must consider, at the very least, the modes of heat transfer that act between the zones of the building, heat gains, material properties, and ambient weather. In addition to this, an EnergyPlus building model will also include components like equipment and occupancy schedules and definitions for individual pieces of equipment. This often equates to thousands of input parameters. Even with highly detailed models, it is difficult to perfectly capture the behavior of a multilevel interconnected system. Hence, BPS is a field that seeks continual improvement that is driven by new research interests and application horizons [13]. 1.2.2 Building Controls On average, approximately 50% of total building energy is used for space conditioning [10]. The operating schedule of these systems is typically determined by setpoints at various stages of the system, where feedback laws govern the supply state of the system. This traditional control method can utilize techniques such as proportional integral derivative feedback to stabilize its operation but they are limited in their ability to achieve optimal performance based on more than one objective. For this reason, older generation buildings have poor energy performance because they were designed primarily to track thermal comfort setpoints. Alternatively, Model Predictive Control (MPC) can be used for complex constrained multivariable control problems and has been implemented extensively in fields such as process control [14], [15], [16], [17]. MPC has become a popular method for controlling a building to achieve optimal comfort standards while also optimizing energy usage, cost, or providing flexibility for demand response. As a general approach to building performance improvement, building control retrofits are a good option because they can reduce operating costs without the necessity of a physical retrofit, provided the necessary data infrastructure is available. Generally, it is becoming increasingly common to have 5 building control and automation systems that provide building managers with the ability to control thermal zones remotely. 1.3 Materials and Resources The code and additional resources that pertain to the content in this thesis can be found online at the GitHub repository listed below. https://github.com/kadenP/TheRoleOfBuildingsInAChangingEnvironmentalEra_PleweMSThesis.git CHAPTER 2 QUANTIFYING EXTERNALITIES IN THE POWER SECTOR The majority of this chapter was extracted from an original publication in the Elsevier SoftwareX Journal, which can be referenced with the following citation: Plewe, K., & Smith, A. D. (2019). PECT: A tool for computing the temporal and spatial variation of externalities related to power generation in the United States. SoftwareX, 9, 61–67. http://doi.org/10.1016/j.softx.2018.12.001 Climate change is expected to have extreme implications within the energy and water infrastructures that we rely on. We already see how these changes are manifested through a changing frequency of extreme events such as droughts, wildfires, and hurricanes [18], [19]. These have the potential to disrupt the systems that supply energy and water to communities and both adaptation and mitigation strategies will need to be adopted to support our current standard of living while sustaining our natural environment. Electricity production accounted for 27.5% of 2017 greenhouse gas emissions in the United States, based on a 2019 EPA inventory of greenhouse gas emissions and sinks [20]. This percentage has been in decline due to the increasing utilization of renewable power sources. Water usage in thermoelectric power plants is also a significant concern in areas with limited water resources. Approximately 87% of electricity generation requires water for cooling in the United States [21], [22]. Accordingly, the power systems and built environment research community need to better understand how energy choices affect the environment through CO2 , NOx and SO2 emissions and water usage [18]. To gain an understanding of the externalities related to our energy systems, both direct and indirect sources need to be considered. Direct emissions and direct water usage occur at the building site, while indirect emissions and water usage often take place far 7 away from the building where central power generation stations are located. Addressing direct emissions involves mainly the transportation sector and on-site combustion while addressing indirect externalities must include the emissions and water use associated with energy conversion, transmission and distribution processes in the power system. Electricity is the largest portion of a typical building’s energy demand in operation, and it is crucial that the associated indirect externalities be quantified. This chapter focuses on methods that can be used to quantify emissions and water usage in the power sector with a high spatial and temporal resolution. Here, high resolution in the spatial sense is anything that can be sourced from individual balancing authorities in regional power transmission systems and high resolution in the temporal domain is anything that is resolved to the hourly or subhourly level. While these highly resolved estimations of emissions and water usage associated with power production aren’t always necessary, it is a knowledge base that is crucial in the operation of emerging technologies such as smart thermostats, smart grids, and automated demand-side management programs. 2.1 Electrical Grid Operation The complexity of electrical generation and distribution systems in the United States necessitates a spatial and temporal modeling approach. Here, fuel mix is defined as the makeup of electric generators that are actively supplying electricity to the grid. Spatial variations in fuel mix have a significant impact on the quantity of emissions emitted, water withdrawn and water consumed for every unit of power supplied to the regional grid. The same can be said for temporal variations in fuel mix, resulting from marginal generators responding to demand and intermittent renewable energy systems, which are increasingly present on the grid. The importance of these variations has been emphasized in studies that look into the spatial and temporal dependence of power system behavior and proposed mitigation strategies [23], [24], [25], [26]. Siler-Evans et al. showed that the benefits of wind or solar vary significantly across regions and that the individual health, environmental and climate benefits do not always coincide when compared in different regions [24]. The importance of temporal and spatial dependency in emission quantification is exemplified in a study on the environmental impacts of electric vehicles [26]. The study cites changes, through time and place, in the fuel 8 mix of the electricity being used for charging as an important indicator of the emissions offsets associated with electric vehicles. Such cases demonstrate the importance of having tools available for researchers and decision-makers that quantify spatial and temporal electricity generation externalities. Methods for capturing changes in regional emission rates through time and space have been previously explored and used in various studies that take advantage of available data sets containing generation data with spatial and temporal resolution [23], [26], [27]. In North America, spatial dependency corresponding to electricity generation is typically attributed to geographical regions that fall within one of four major interconnections: Eastern Interconnection, Western Interconnection, Texas Interconnection, and Quebec Interconnection. These interconnections and the divided regional entities within are shown in Figure 2.1. The management of these interconnections and regional entities is managed by the North American Electric Reliability Corporation (NERC) and carried out by independent system operators (ISO), regional transmission organizations (RTO), and/or balancing authorities. Figure 2.1. Map of NERC interconnections and regional entities [4] 9 The fuel mix can be quantified spatially by mapping independent system operators (ISO), regional transmission organizations (RTO), and/or balancing authorities. ISO and RTO entities are independent organizations that facilitate the trade of electricity as well as control and monitor activity on regional electrical grids, and balancing authorities have the specific task of managing the interplay between supply and demand on the electrical grid [28]. The EPA annually publishes averaged emission and generation data for subregions and balancing authorities in the eGRID database [3]. This database is used by researchers to calculate average emission factors and marginal emission factors–emission factors from marginal providers [24], [26], [29], [30], [31]. It has also been used to calculate region-specific emission factors for generation fuel types [27]. The subregions of the U.S. represented in the eGRID database are shown in Figure 2.2. The tool presented here, PECT, is only useful for studies within U.S. regions that are included in the associated databases. A similar tool could be developed for other regions that have a similar set of publicly available databases. Figure 2.2. Map of the subregions represented in the latest eGRID database published in 2014. Four subregions are excluded in this image: AKMS and AKGD, in Alaska and HIMS and HIOA, in the Hawaiian Islands. Crosshatching indicates that the geographical area covers multiple eGRID regions [3]. 10 In the United States, the EPA collects Continuous Emission Monitoring System (CEMS) data in hourly intervals for power generating plants with capacities greater than 25 MW. The data can be accessed through an online platform provided by the EPA [32], and have been used in parallel with the eGRID database to study spatial and temporal emission trends in the U.S. [33], [34]. A study done at Carnegie Mellon University used CEMS data at the facility level from the EPA’s open web portal [32] and state level generation and fuel consumption data from the EIA’s open web portal [35] to estimate carbon intensity for annual, quarterly and monthly time-frames with different spatial resolutions [36]. Shivley provides a method of emission quantification that is easily reproducible and open source, much like the work being presented here [36]. The time dependency of emissions coming from power generation can be attributed to changes in marginal load providers such as generating plants that operate to meet marginal demand and intermittent suppliers such as solar and wind. Past studies have utilized unit commitment and dispatch data [23], locational marginal pricing data [25], [27], explicit fuel mix data [25], [26] and explicit emissions data [24], [29], [30], [33], [34], [36] to capture the time dependence of emission factors. The same spatial and temporal dependency can be applied to water usage in power generation systems. Power generation water usage is quantified in two categories: water consumption and water withdrawal. Water consumption is water that is permanently removed from a natural water source or utility, typically attributed to closed-loop cooling systems in thermoelectric power plants or evaporation on reservoirs behind hydroelectric dams. Water withdrawal describes all water that is removed from a source, often attributed to once-through cooling systems in thermoelectric power plants [37], [8]. Forty-five percent of water withdrawals for all uses was attributed to thermoelectric power generation in 2010 [7]. In 2011, thermoelectric plants with once-through cooling systems supplied approximately 23% of electricity in the United States and accounted for 64% of the water withdrawn for power generation. In the same year, thermoelectric plants with closed-loop cooling systems accounted for 35% of electricity generation and 17% and 88% of water withdrawn and consumed, respectively, by electricity providers in the United States [8]. Water consumption and withdrawal in the power sector have received scientific attention in an effort to understand the water-energy nexus and its future implications. Multiple 11 studies have provided long term case scenarios that give insight into future expectations for the interplay between water and energy systems under various scenarios and in regions with various levels of water scarcity [38], [39], [40], [41]. Peer et al. present a method for obtaining spatial and temporally resolved water consumption data using water consumption rates [23]. Peck and Smith model electricity generation and water usage on the plant level to quantify the water consumption and withdrawal and water resource impact under different generation scenarios [21]. All of these studies cite the importance of avoiding future energy development that inherently neglects to consider the water-energy interplay on a locational basis. Such considerations will be integral in the longevity of our energy infrastructure as we experience climatic change [8], [22], [38], [41], [42], [43], [44]. 2.2 Power Externality Correlation Tool The acquisition of temporally resolved fuel-mix and environmental externality data for various locations in the United States currently requires the combination of multiple data sources. Tools that provide data that can be used in environmental analysis are often location limited [27], have low temporal resolution or only provide historical data [3], [36], [45], [46]. These limitations make it difficult for environmental analyses for power purchases to be performed on a cross-disciplinary platform. In an effort to address the limitations in the available tools and resources, the Power Externality Correlation Tool (PECT) allows researchers and policy makers to obtain the following data either in real-time or historically for various regions within the United States: • Fuel mix at hourly intervals. • Total CO2 , NO2 and SOx emissions at hourly intervals. • Water consumption and water withdrawal at hourly intervals. • Fuel-specific CO2 , NO2 , SOx , water consumption and water withdrawal factors. 2.2.1 Materials and Methods PECT takes advantage of external databases and application programming interfaces (API) for geographical positioning, electrical generation fuel mix data, water factors and emission factors. The WattTime API [1] was used to obtain locational electricity generation 12 data for different fuel types (coal, natural gas, nuclear, biogas, wind, geothermal, solar, hydro and biomass). WattTime collects generation data for balancing authorities and power markets that supply a consistent stream of data. The balancing authorities and markets that are available through the software are listed in Table 2.1. RT5M and RTHR markets contain historical data at five-minute and one-hour intervals, respectively. The DAHR market contains day ahead forecasted data at one-hour intervals. Electricity generation information from WattTime is received in a JSON dictionary format with the time stamp, total carbon output, generation by fuel type, balancing authority and market. Data are provided at hourly and five-minute intervals, and here it is averaged and reduced to hourly intervals. Version two of the WattTime API was released in 2018 and includes a larger selection of spatial data across the U.S. and Europe, but PECT has not yet been updated to take advantage of this functionality. Water usage estimates are obtained using water factors that have been compiled in a comprehensive literature review [2]. The accuracy of these factors will decrease as time goes on due to changes in the technology that is used in operable power plants. The factors used in this work are presented in Table 2.2. One assumption worth noting that pertains to Table 2.2 is that all withdrawal and consumption factors, aside from geothermal, are for plants that utilize cooling towers. Since the available data for geothermal plants were limited [2], the geothermal withdrawal and consumption factors were taken as an average of all the technologies listed. Water consumption and withdrawal are calculated using Table 2.1. Generation Data Availability from WattTime API [1] Balancing Authority BPA CAISO CAISO CAISO ERCOT ISONE MISO MISO PJM PJM SPP Market RT5M DAHR RT5M RTHR RTHR RT5M DAHR RT5M RT5M RTHR RT5M Available since (UTC) Feb. 26, 2014 Aug. 5, 2014 Feb. 26, 2014 Feb. 1, 2013 Feb. 26, 2014 Dec. 29, 2013 Dec. 2, 2015 Feb. 26, 2014 Feb. 27, 2014 Sept. 24, 2016 Feb. 25, 2014 13 Table 2.2. Water Consumption and Withdrawal Factors [2] Fuel Type Coal Natural Gas Nuclear Biogas Wind Geothermal Solar Thermal Solar PV Hydro Biomass Water Consumption Factor [gal/MWh] 687 205 672 235 0 251 786 0 4491 235 Water Withdrawal Factor [gal/MWh] 1005 225 1101 878 0 251 786 0 0 878 10 WC = ∑ En × WFc,n (2.1) n =1 10 WW = ∑ En × WFw,n (2.2) n =1 where WC is the water consumed, WW is the water withdrawn, En is the electrical energy produced using fuel type n (10 fuel types are considered here), and WFc,n and WFw,n represent the water consumption factor and water withdrawal factor, respectively, for fuel type n. Emission estimates for CO2 , SO2 and NOx are included in PECT. CO2 , SO2 and NOx emission factors were obtained from the EPA’s eGRID database [3] for each fuel type and are unique to each balancing authority included in the model. To obtain these fuel-based emission factors, a weighted average of the emission factors for plants of the same fuel type is taken for the balancing authority determined from the user-input location. This is done using EFn = ∑im=1 EFi Gi ∑im=1 Gi (2.3) where EFn is the emission factor for fuel type n, EFi is the emission factor for plant i and Gi is the total annual electrical energy generated for plant i located in a balancing authority region. The emission factors for the regions supported by the WattTime API are listed in 14 Table 2.3. In this table and subsequent calculations, gas emission factors are used for both natural gas and biogas. 2.2.2 Software Architecture The basic architecture for PECT is depicted in Figure 2.3. PECT was developed using Python 3.6 and is provided as an open-source script at the Site-Specific Energy Systems Laboratory GitHub repository [47]. PECT provides a simple means of accessing data on the Table 2.3. Calculated emission Factors for CO2 , NOx and SO2 using Equation 2.3. Emission factors are calculated with eGRID data [3] for Coal, Gas and Biomass for each balancing authority included in the WattTime API [1] Fuel Type BPA Coal Gas Biomass CAISO Coal Gas Biomass ERCOT Coal Gas Biomass ISONE Coal Gas Biomass MISO Coal Gas Biomass PJM Coal Gas Biomass SPP Coal Gas Biomass CO2 Emission Factor [lbs/MWh] NOx Emission Factor [lbs/MWh] SO2 Emission Factor [lbs/MWh] 2327 925.9 84.13 2.201 0.3362 2.684 2.119 0.004868 1.244 2147 846.6 117.1 0.7201 0.2016 2.962 3.701 0.005998 0.3640 2238 864.5 – 1.145 0.3421 – 4.749 0.005348 0.3064 2075 886.4 801.6 1.691 0.2092 3.239 2.000 0.05643 0.8604 2196 931.5 296.3 1.816 0.6506 3.545 4.429 0.07093 1.674 2088 934.6 733.2 1.998 0.2348 4.386 5.017 0.02267 1.431 2212 1032 145.5 1.798 0.9202 2.094 3.654 0.03365 1.735 15 Figure 2.3. PECT architecture outlining the process that is used to obtain emissions and water usage estimates. 16 WattTime API Database and calculating additional electrical generation metrics. Required user inputs include desired location, time frame and WattTime API credentials. To use PECT, the user must have a WattTime API account, which can be obtained for free at the WattTime website [1]. The Google Geocoding API is used to convert the desired location to geographical coordinates. These coordinates and the specified time frame are used for the WattTime query to download the generation fuel mix data. With the generation data in a JSON dictionary format, the relevant data points are extracted and water consumption, water withdrawal, CO2 emissions, NOx emissions and SO2 emissions are calculated with their associated factors. Results are provided to the user in a csv spreadsheet format, which can then be used for graphical presentation or additional analysis in a number of applications. 2.2.3 Software Functionalities The main function of PECT is to provide emissions and water usage estimates. Users with varying levels of experience, including those with no experience developing in Python, can take advantage of PECT either by downloading the materials and using it in its current state or by using it as a building block for separate Python applications. It can also be operated in Google Collaboratory on a virtual kernal, eliminating the need for the user to download software and dependencies. In addition to externality quantification, PECT can be used to obtain emission and water factors related to specific fuel types and locations. PECT is dependent on the availability and validity of the data reported by ISOs. Table 2.1 lists the entities that report to WattTime. However, this does not guarantee that all electricity generators are included or that the data sets are free from temporal discontinuities. For the examples included in this work, a simple interpolation scheme was used to fill gaps in data. The numbers used here were checked through comparisson with well known databases such as the eGRID database [3], Power Profiler [45], and the EIA electricity data browser [46]. 2.3 Case Studies PECT writes the results to a csv output file. The output data include generation data for each fuel type, emissions data, water consumption data, water withdrawal data, 17 emission factors, water consumption and withdrawal factors. A sample output is included in the PECT GitHub repository. From this data set, studies can be performed and the environmental impacts of electricity purchases can be represented visually. As an example, Figure 2.4 shows the percentage of generation, emissions and water usage attributed to each fuel type reported for a location in the MISO region. This was obtained directly from the output data using Pn,t = Xn,t × 100 Xt (2.4) where Pn,t is the percentage of a metric (generation, emissions, etc.) for fuel type n at time t and Xn,t and Xt are the metrics for fuel type n and total metric, respectively, at time t. To demonstrate the regional resolution supplied by PECT, the generation fuel mix obtained for three different locations is shown in Figure 2.5. The MISO, PJM, and ISONE Figure 2.4. Total generation, emissions and water usage attributed to various fuel types reported by the MISO balancing authority. Software query was performed for Kalamazoo, Michigan for April 22, 2018. The reg. fit lines are obtained using relaxed least squares regression. 18 Figure 2.5. Generation mix for three locations: (a) Michigan (top), (b) Virginia (bottom left), and (c) Vermont (bottom right) on July 17, 2017. The balancing authorities for these locations are (a) MISO, (b) PJM, and (c) ISONE. regions were chosen for this example because they were found to be well-populated with fuel mix data during manual testing of the API calls in PECT. The locations were selected by passing in the names of cities within each region (i.e., Kalamazoo for MISO). Many of the other regions in the United States either do not report fuel mix data or ignore major energy producers, making the associated fuel mix data less accurate. The implications of purchasing electricity in different regions at the same time can be explained using Figure 2.5. Referencing the fuel mix that is shown for Michigan, Virginia and Vermont for a single day on July 17, 2017, it is clear that the emissions associated with a single unit of electricity purchased at 12:00, for example, would be the greatest in Michigan and smallest in Vermont. This is because electricity generation in Michigan has a larger percentage of generation from coal fired power plants compared to both Virginia and Vermont using more nuclear and natural gas at the same point in time. From an energy management standpoint, assuming that emissions and water usage are a mitigating factor, this information could influence energy management decisions in various regions. For example, the environmental considerations of electric versus natural gas heating in 19 buildings would vary significantly in Vermont and Michigan. The implications of purchasing electricity at different times in the same region can be explained using Figure 2.4 and Figure 2.5. Both figures show the temporal variations in the fuel mix within individual regions. Looking at the plot of Vermont in Figure 2.5 as an example, the percentage of electricity coming from nuclear and natural gas changes by approximately 10% throughout the day. This information is particularly useful in cases that employ demand shifting strategies with the goal of reducing emissions. For example, shifting loads from the end of the day to the beginning of the day on July 17, 2017 in Vermont would reduce emissions because the percentage of electricity coming from nuclear is highest between the hours of 1:00 A.M. to 12:00 P.M. 2.4 Implications PECT contributes to the power systems and building energy systems research communities, policy development personnel, and conscientious energy managers by providing a simple method for quantifying spatial and temporal externalities that come from energy usage. Previous work has focused on developing methods to quantify energy system externalities, but tends to focus on individual regions and does not provide comprehensive tools to the general public [23]. PECT extends existing methods for quantifying these externalities and integrates them into a tool that is applicable to multiple regions in the United States with the purpose of supporting interdisciplinary efforts to better manage our power systems. While water usage and emissions have been modeled separately by many research groups, the study of both together has been limited. By providing a means to quantify water usage and emissions side-by-side, PECT enables a direct comparison of the environmental costs and benefits of various fuel mix scenarios. Comparing water usage and emissions side-by-side allows researchers and policy makers to avoid implementing strategies that may have unintended consequences because of conflicting environmental values, as it has been shown that water conservation and pollution reduction do not always coincide [23]. Carbon accounting is important for decision making on many societal platforms: governmental, industrial, commercial, residential, etc. Just as carbon emissions are closely monitored in the United States and across the globe, other polluting emissions and water 20 usage can be monitored to inform the public, researchers and policy makers so that they are better informed when making benchmarks for future environmental goals [48]. This is especially important in areas that expect to significantly increase generating capacity to meet growing demand because it can influence the design of generation plants (fuel type, cooling system, placement, etc.). PECT is an important example of how researchers can synthesize and utilize available data to unveil valuable information for policy makers and the general public. The open source development of tools such as this one supports future development of more complicated models that could provide a sound basis for data-driven environmentally-minded decision making in many different sectors. 2.5 Summary Models that quantify the externalities related to power systems have received a high level of interest in recent years because of the implications that these externalities have within our communities and the environment. This paper briefly reviewed some of the literature that has been published that introduces and builds upon such models. A common goal within this research community is to identify methods that best represent externalities with spatial and temporal resolution. This is important because power systems are diverse in both time and place across the United States and other parts of the developed world. The models introduced in literature describe the methods used in detail but rarely focus on developing tools that implement these methods for general use. Such tools are necessary because decision makers often lack the time and resources to explore and implement available methods. PECT is an important step in implementing practical methods to calculate metrics that provide insight into the environmental consequences of power generation. The additional consideration of power prices and a full life-cycle analysis of power generation would increase the usefulness of PECT. In particular, showing cost metrics side-by-side with fuel mix, emissions and water usage data would significantly increase the applicability of PECT within the U.S. economy. CHAPTER 3 MODEL PREDICTIVE CONTROL IN BUILDING THERMAL SYSTEMS The energy systems that are used to maintain thermal comfort and air quality within a building have always been a crucial part of modern life. They play a major role in maintaining productivity, comfort, and health and are now regulated and standardized in many countries [49]. With the increased adoption of heating, ventilation and air conditioning (HVAC) and other energy systems in buildings around the world, the buildings sector now makes up as much as 40% of total source energy consumption in developed countries [9], [10]. A further breakdown of this energy usage shows that HVAC systems typically account for more than 50% of a building’s total energy demand [10]. This indicates that building HVAC systems should be a major focus in efforts that aim to mitigate the impacts associated with global energy demand. Improvement of building HVAC controls is a prime option for a building energy reduction strategy because there typically isn’t a need for large renovations. They can also present a multitude of advantages as building systems become more integrated at the single building scale and grid-scale due to their ability to accomplish load shifting, peak shaving and demand reduction without compromising occupant thermal comfort. Most buildings are now equipped with some version of a Building Automation System (BAS) which provides the capability for remote access and control of various components of building energy systems. These BAS’s can automate the control of building systems based on a strict set of rules and logic. With this capability and the improved access to operational data, building energy managers now have the option to adopt more sophisticated control strategies with a relatively low implementation burden with regard to physical infrastructure [50]. Model-based control strategies have grown in popularity in the buildings sector in parallel with advancements in building modeling [51], [52]. These allow for improved 22 control decisions that are optimized based on a preordained knowledge of how they will affect different aspects of building operation. In addition to the advancement of building energy simulation techniques, improved weather forecasting has increased the accuracy of building energy forecasting since this is one of the most influential disturbance variables in building energy systems [53], [54], [55], [56]. Model predictive control (MPC) is a control methodology that takes advantage of these two things, a building model and disturbance forecasts, to optimize building control based on how the building is expected to behave over a given time horizon. The performance of an MPC scheme is dependent on the accuracy of the model and disturbance predictions, which typically consist of occupancy schedules and local weather conditions [57], [53]. Hence, it is important to acknowledge the difference between results obtained from simulated applications and real building case studies. MPC still faces many challenges with regard to real building implementation, which mainly pertain to the complexity of accurately modeling a building for real-time control [58], [59]. The application of building modeling in real time control is currently at the forefront of building energy simulation research [13] which will significantly increase the rate of adoption of MPC in real building control. 3.1 Problem Definition In this study, an MPC technique is presented in a simulated small office building that is modeled based on the Department of Energy’s small office prototype building [60]. For this application, the building model is reduced to a simple thermal resistance-capacitance network with heating capabilities. This method for modeling building energy devices and thermal transfer has been utilized extensively for similar control applications [61], [62], [63], [64], [65]. Modeling the building in this way allows the building model to be used in a simple convex optimization problem in the MPC control algorithm. The traditional method for applying MPC in building control utilizes a medium-term prediction horizon (greater than 24 hours) because this allows the passive thermal properties of the building to be utilized in the optimal control when accounting for large ambient temperature fluctuations. However, since the building model used in this case neglects some of the most important thermal storage components in the building, passive thermal design is not a focus here. In this chapter, a relatively short forecasting horizon is first used to test the performance 23 of an MPC algorithm based on energy usage and thermal comfort without the added benefits that come with large forecast horizons. Three different control techniques are presented for comparison: a basic on/off controller, an MPC controller with on/off heat input and an MPC controller with variable heat input. These short forecast horizon simulations use a time step of one minute and don’t consider variations in occupancy within the building. As a follow-up, the variable heating MPC technique is used in a case where a 15-minute time step and 24-hour forecast horizon is used in a simulation that considers the occupancy schedule defined by the DOE in the small office reference building. One of the benefits of MPC is demonstrated in this case, as it is used as a method of peak shaving for high demand intervals. 3.2 System Model The system that is being modeled in this study is the small office reference building that is part of the DOE commercial prototype building model database [60]. The building envelope and interior are shown in Figure 3.1. The building has five conditioned thermal zones that interact with each other through internal walls, four of which are external and exchange heat with the surrounding environment through external walls. In addition to heat exchange with spaces horizontally, each thermal zone exchanges heat with the attic of the building. A simplified thermal RC model is created based on the layout and geometry of this building and does not perfectly reflect the EnergyPlus building model due to major simplifications that were used in order to reduce the analytical complexity of the state space model. While this study is not prioritizing a perfect representation of the building physics, it is important to acknowledge some of the major assumptions that were made in the creation of the thermal RC network model. These are listed below. • The floor acts as an adiabatic boundary. • The thermal capacitance of the walls and all other building materials is neglected. • The walls are modeled as a homogeneous medium with constant thermal resistance. • Only heat transfer due to conduction between spaces is considered. 24 Figure 3.1. Isometric and top view of the interior of the small office prototype building used in this study. The outdoor temperature data that are used for this analysis were extracted from a TMY3 weather data file for Lansing, MI., which is in IECC climate zone 5A [66]. The temperature was obtained at one-minute intervals for an entire year, and simulations were performed over winter periods in January. Over this period, the building will only demand heating, which is more practical for the thermal RC circuit model since air conditioning provides additional complexity that would otherwise be neglected here. The design heating loads that are calculated in the EnergyPlus simulation were used here for scaling purposes when selecting the heat input power for on/off control schemes. These design load calculations are provided in Table 3.1. 3.2.1 Thermal RC Model The small office building was modeled here using a single thermal resistance for the barriers separating spaces and a single capacitance for the air contained within a space. This is considered a 1R1C network and is an overly simplified model because of the assumptions mentioned earlier. The circuits that represent the connections between each of the thermal Table 3.1. Heating design power calculated by EnergyPlus for each thermally regulated zone for a design temperature of -16.30 ◦ C. Zone 1 3323.29 EnergyPlus Design Power (W) Zone 2 Zone 3 Zone 4 1940.61 2942.91 1925.84 Zone 5 1325.36 25 zones considered here, the external environment and the control heat input are shown in Figure 3.2. We start the formulation of the state space equation for this thermal model with an energy balance for each zone which is written as ρVi C Ti (k + 1) − Ti (k ) = ∆t p ∑ j =1 Tj (k ) − Ti (k ) + qi ( k ) Rij (3.1) where i = 1, 2, ..., 6 is for each thermal zone (including the attic which is denoted here by subscript a or the number 6) which has a unique volume, Vi , and heat capacity, C. In this model, however, the total capacitance for some of the zones is identical due to symmetry. A finite difference method is used here to represent the change in temperature, T, within a zone over a time step, ∆t. On the right-hand side of Equation 3.1, the summation accounts for heat transfer between thermal zone i and all of the adjacent spaces j = 1, 2, ..., p including outside for all of the zones except zone 5. The total thermal resistance, Rij , accounts for the area and thermal resistance of the barrier separating the two adjacent zones and has the units of K/W. Finally, qi (k ) is the heat power input to the thermal zone at time step k. (a) External Zone Circuit (b) Core Zone Circuit (c) Attic Circuit Figure 3.2. Thermal RC circuits for external and core zones with heating input and the attic zone. 26 After doing this heat balance for each thermal zone, a state space representation for the buildings thermal network is constructed such that C T ( k + 1) − T ( k ) = AT (k) + Bq(k) + DTo (k) ∆t (3.2) where C ∈ R6×6 , A ∈ R6×6 , B ∈ R6×6 , and D ∈ R6 now contain the information that completely defines the dynamic system with zone temperatures T ∈ R6 . On the left hand side of Equation 3.2, C contains the total thermal capacitance for each of the separate zones in units of J/K and is calculated as  V1   0 C = ρC   ..  .  0 0 ··· 0    0  ..  .  · · · Va V2 · · · .. . . . . 0 (3.3) On the right hand side of Equation 3.2, A contains the information relating each of the building zone temperatures contained in T as the inverse of the connecting resistances with units of W/K, B indicates which zones are thermally regulated with a heat input q ∈ R6 and D relates the building zone temperatures to the outdoor temperature To (k ) ∈ R1 . The formulation of A, B, and D is shown below.  1 1 0 R1j R12 − j∑ 6 = 3   1 1  − ∑ R12j R23  R21 j 6 =4   1  − ∑ R13j 0  R32  j 6 =1 A=  1 1 0  R41 R34    1 1  R151 R52 R53    1 1 1 R1a R2a R3a    1 0 ··· 0      0 1 · · · 0   B = . . . ..  , D = . . . . . . .   0 0 ··· 0 1  R1o   1   R2o     1  R   3o   1  R   4o     0    1 R ao 1 R14 1 R15 1 R1a 0 1 R25 1 R2a 1 R34 1 R35 1 R3a 1 R45 1 R4a 1 R54 − ∑ R15j 1 R5a 1 R4a 1 R5a − ∑ R11j −∑ j 6 =2 1 R4j j j            ,          (3.4) 27 where the subscripts a and o denote the attic temperature and the outdoor temperature, respectively, and the numeric subscripts denote the zone number as listed in Figure 3.1. The ij pair represents the connection between zone i and j and Rij = R ji , so A is a symetric matrix. In Equation 3.2, matrices C, B and D can be expanded for an arbitrary number of zones trivially by considering the volume, heat input, and its adjacency to outside for each zone. The A matrix can be generalized for an arbitrary number of zones, M, by first defining sets, {S1 , S2 , ..., S M }, that contain the zone numbers for every zone that they exchange heat with (i.e., adjacent zones). The A matrix can then be defined as Aij =  − ∑    q ∈ Si    1 Riq i=j 1 Rij i 6= j and j ∈ Si 0 else (3.5) for any size, A ∈ R M× M , where i and j are the row and column indices, respectively, for A. Further reducing this in order to obtain an explicit expression for the state response variables, x (k + 1), in terms of the previous states, x (k ), the inputs, u(k ), and the disturbance variable, w(k ), we get x (k + 1) = Âx (k ) + B̂u(k ) + D̂w(k) (3.6) Â = I + ∆tC−1 A, B̂ = ∆tC−1 B, D̂ = ∆tC−1 D (3.7) where and I is the identity matrix, I ∈ R6×6 . This is the final form of the state space equation that is used to relate the temperature of each of the thermal zones to the heat input and the outdoor temperature disturbance. All of the resistance and capacitance values were obtained using data from the EnergyPlus building model to preserve the relevance of the results obtained here within the actual system. The capacitance and resistance values that are used are listed in Table 3.2 and Table 3.3, respectively. To verify the dynamics of this model, an isolated simulation was performed where the heat input, u, was set to zero and the initial temperature of all of the thermal zones was set to the initial ambient temperature. With these initial conditions, we expect the interior temperatures to track the ambient temperature with a temporal lag that is dependent on 28 Table 3.2. Data used to calculate the thermal capacitance of each zone. ρ (m3 /kg) 1.205 C (J/kgK) 1005 V1 (m3 ) 346.5 V2 (m3 ) 205.3 V3 (m3 ) 346.0 V4 (m3 ) 205.3 V5 (m3 ) 456.5 Va (m3 ) 383.4 Table 3.3. Total thermal resistance values obtained using resistance and areas listed in the EnergyPlus small office model. Values are oriented to represent the resistance between the row and column nodes. Zone 1 Zone 2 Zone 3 Zone 4 Zone 5 Attic Zone 2 0.04547 - Total Thermal Resistance (K/W) Zone 3 Zone 4 Zone 5 0.04547 0.01817 0.04547 0.03800 0.04547 0.01817 0.03800 - Attic 0.00864 0.01457 0.00864 0.01457 0.00655 - Outside 0.03262 0.04893 0.03262 0.04893 0.00477 the thermal capacitance of each zone. The result of this simulation is shown in Figure 3.3. This behavior aligns with the physical expectation, which verifies the model formulation shown here. 3.3 Methodology In this problem, MPC is applied as a constrained tracking problem using the thermal building model introduced in Section 3.2.1. The system is considered stable, as the physical relationship between the system and the environment necessitates that the state of the building thermal zones follow the outdoor temperature as time goes to infinity. This is the case over the time horizon that is being considered here based on the demonstration shown in Figure 3.3. Three different methods for controlling the temperature to a setpoint of 22 ◦ C were tested. The first is a basic on/off control scheme where a single rated heat input is applied to a zone whenever it drops below the setpoint temperature. It is a live controller that has no forecasting ability and only uses binary feedback to set the state of the heat input. The second is an MPC scheme that utilizes a single rated heat input, functioning similar to the basic on/off scheme with forecasting and optimization capabilities. The third is an MPC scheme that utilizes a variable heat supply system, allowing the MPC optimization to 29 Figure 3.3. Outdoor dry bulb temperature for a 48-hour period in January with the thermal response of the zone temperatures without heat input. directly control the magnitude of the heat being applied to each thermally regulated zone at every time step. Two different preprogrammed setpoint schedules are considered. One is constant, representing a case where the temperature of the office building requires thermal conditioning 24 hours a day. The other is a variable temperature band schedule that is the same as the one designed for nonholiday weekdays in the EnergyPlus small office model. This schedule represents a typical workday occupancy schedule. The temperature setpoints and heating and cooling bands used are shown in Table 3.4. Table 3.4. Heating and cooling setpoint schedules used over the winter control simulation period. Hour 0-6 6-7 7-8 8-18 18-19 19-24 Constant Setpoint Simulations Tsp (◦ C) 22 22 22 22 22 22 Variable Temperature Band Simulations Heating Setpoint (◦ C) Cooling Setpoint (◦ C) 15.56 29.44 17.8 27.8 20 25.6 21.11 23.89 21.11 29.44 15.56 29.44 30 3.3.1 Basic On/Off Control The first method was a simple on/off scheme where the heat input to each thermally regulated zone could either be zero or some rated power. The constant heating power supply was chosen to meet the design loads listed in Table 3.1. Because Zone 1 and Zone 3 are identical in the simplified thermal model, the same constant heating power was used for these zones. This is also the case for Zone 2 and Zone 4. The constant heating input values that are used for this control scheme are listed in Table 3.5. This type of on/off control is the most basic form of thermostat temperature regulation. It would not be fully practical as it is implemented here because it neglects time constraints for ramp-up periods and has no method for facilitating steady operation. For a better comparison between MPC and standard control techniques, more effort would need to be put into developing the baseline controller. For example, on/off changes could be controlled based on a temperature band instead of a single setpoint as is typically done in real buildings. Regardless, the algorithm used to carry out this basic on/off control technique is shown in Table 3.6. It controls the heat input based only on the placement of the state variables relative to the setpoints. It doesn’t consider the steady state error or error magnitude, and as a result, the heating supply is relatively unstable and switches from on to off at a high frequency. This would be undesirable in a real building scenario because 1) thermal comfort can be reduced as a result of these fluctuations, 2) the startup times of some heating systems limit the frequency at which on/off cycles can occur, and 3) the lifetime and efficiency of a heating system is affected by the on/off cycle frequency. Because these things aren’t accounted for quantitatively in the results presented here, the performance of the on/off control method is overly optimistic. Table 3.5. Constant heat supply levels that are used for on/off control schemes. Zone 1 4000 Heat Supply (W) Zone 2 Zone 3 Zone 4 3000 4000 3000 Zone 5 2500 31 Table 3.6. Basic on/off control algorithm for each time step, k, in simulation do if zone temperature is less than the setpoint then apply heating power to zone else set heating power to zero end calculate zone temperatures at next time step end 3.3.2 Model Predictive Control Two versions of MPC were tested for this problem, one utilizing a single rated heat supply for each conditioned zone as was used in the basic on/off control technique and one that utilizes a variable heat supply. Both methods use the same optimization algorithm with the only difference being the constraint on the heat supply. Within the context of a building energy system, the main goal of a control system is to maintain thermal comfort. Beyond this, it is desirable to minimize the associated cost of running the heating system which can be evaluated by energy usage, emissions, or monetary cost, for example. With this basis, the MPC optimization used here considers the heat input sequence, the difference between the zone temperature and the setpoint, and rate of change of the heat input. Indirect emissions are not explicitely accounted for here, but they are influenced by demand and total energy usage. The following general problem formulation is proposed for this use case: min u(k|k ),u(k +i |k),...,u(k + N |k ) J1 ( u ) + J2 ( u ) + J3 ( u ) s.t. 0 ≤ u ≤ u x≤x≤x (3.8) where u ∈ R6 is the heat input for each zone and x ∈ R6 is the temperature for each zone. J1 (u) is the cost associated with the magnitude of the heat input to each thermally regulated zone and J2 (u) is the cost associated with the deviation of the zone temperatures from the setpoint temperature. The third cost, J3 (u), accounts for the difference in heat input at two sequential time steps and is utilized to reduce sharp demand peaks. The complete expressions for J1 (u), J2 (u) and J3 (u) are provided later in Equations 3.11, 3.12 32 and 3.13, respectively. All of the objective functions are represented only as a function of the input, u, that is being manipulated in the optimization procedure. However, J2 (u) is also inherently dependent on the disturbance temperature, w, and the state temperatures, x. Since we are posing this problem in a season that only utilizes heating, the input constraints are set to maintain the heat input between zero and some maximum value. In addition to applying a cost to zone temperature deviation from the setpoint, a hard constraint for the state variables is applied to ensure that the temperature is maintained within a comfortable range. The forecasting horizon is defined by N and sets the prediction horizon for the zone temperatures based on the known outdoor temperature. In order to consider the predicted temperatures in the objective function and constraints, an augmented state space representation is formulated to relate the temperature at a future time step, k + i, for i = 1, 2, ..., N to the temperature at the current time step, k. This augmented state space representation is defined as x (k + i |k) = M x (k|k) + C u(k + i |k) + D w(k + i ) (3.9) where the augmented state matrices, M, C and D are constructed by expanding the original state space form shown in Equation 3.2 for k = 1, 2, ..., N.    B̂ 0 Â  Â2   Â B̂ B̂    M =  . , C =  . ..  ..   .. . Â N −1 B̂ Â N −2 B̂ Â N  D̂ 0  Â D̂ D̂  D= . .. ..  . Â N −1 D̂ Â N −2 D̂  ··· 0 · · · 0  ., .. . ..  · · · B̂  ··· 0 ··· 0  ..  .. . . · · · D̂ (3.10) With this augmented state space representation, we can operate an MPC algorithm that obtains the optimal heat input for each time step in our simulation, [û∗ (1), û∗ (2), ..., û∗ (k end )] T . The algorithm that is used here for a simulation over a time sequence, k = 1, 2, ..., k end consists of two major steps. 1) The optimal control sequence over a forecast horizon, N, is obtained by evaluating the formulated cost functions and constraints, and 2) the optimal control sequence for the current time step is extracted and applied. This process is then stepped forward and repeated until k end is reached. Here, the control laws are applied to the 33 building model. This is different from how MPC would be implemented in a real building, where the optimal control laws are applied directly to the real building system at each time step and the actual state response values for the current time step is read from sensors as feedback. A more detailed description of this algorithm is depicted in Table 3.7. The operation of the MPC control algorithm is determined largely by the formulation of the cost function and constraints that are used in the optimization step. The cost function that was formulated here was split into three different parts. The first weights the importance of minimizing the heat input, J1 (u), and we refer to it as the input cost. The second part weights the importance of minimizing the deviation of the zone temperature from the setpoint, J2 (u), and we refer to it here as the thermal comfort cost. The third is referred to as a differential cost and captures the importance of reducing sharp demand peaks that can occur if the system attempts to correct state errors over a short interval by calculating the difference between sequential heat input signals. The importance of this can vary based on the type of heating system being used or the demand response responsibilities of the building energy manager. Since this is something that varies in importance, it is one of the main focuses of the simulations presented in Section 3.4.2. The mathematical expressions for J1 (u), J2 (u) and J3 (u) are Table 3.7. Model Predictive Control algorithm used for on/off and variable heating supply MPC simulations. for each time step, k, in simulation do calculate the optimal control sequence over forecast horizon (N), {u∗ (k|k), u∗ (k + 1|k), ..., u∗ (k + N |k)} if variable heating control is available then apply optimal heat input for the current time step, û∗ (k ) = u∗ (k |k ) else if heat input at k greater than zero, u∗ (k) > 0 then turn heat on, û(k ) = qconst else turn heat off, û(k ) = 0 end end calculate the zone temperatures for next time step using Equation 3.2, x (k + 1) = Âx (k ) + B̂u(k ) + D̂w(k ) end 34 J1 (u) = uT Qu (3.11) J2 (u) = (x − Tsp )T R(x − Tsp ) (3.12) J3 (u) = (u − D̃u)T Q̃(u − D̃u) (3.13) where Q ∈ R6N ×6N , R ∈ R6N ×6N and Q̃ ∈ R6N ×6N are weighting matrices that determine the importance of each objective and D̃ ∈ R6N ×6N is a matrix that facilitates the calculation of a temporal difference in heat inputs for each zone. Q = ρI, R = ψI, Q̃ = φI  0 0   .. .  D̃ =  0 0   .. . 0 0 0 0 0 0 .. .. .. . . . 0 0 0 0 0 0 .. .. .. . . . 0 0 1 0 0 0 .. .. .. . . . 0 0 0 0 0 0 .. .. .. . . . 0 ··· 1 ··· .. . . . . 0 ··· 0 ··· .. . . . . (3.14)  0 0  ..  .  1  0  ..  . (3.15) 0 0 0 0 0 0 0 0 ··· 0 where I ∈ R6N ×6N is the identity matrix and ρ, ψ and φ are constants that can be varied to change the weight of each cost. Here, ρ = 0.001, ψ = 5000 and φ = 0 is used for the first set of simulations. In the second set of simulations, a parametric analysis is performed by varying φ with the same values of ρ and ψ. These parameters act as tuning devices that allow the designer of the controller to set the importance of each objective for specific use cases. The weight of the objectives at individual time steps over a single forecast horizon could be varied by defining the weights ρ, ψ and φ as vectors instead of setting them as constant values. Assigning different weights to the objectives at individual time steps could be useful because occupancy and weather predictions become more uncertain as you look further into the future. So, for example, it would make sense to scale the weight of objectives that depend on occupancy or weather forecasts based on temporal proximity. The objective functions consider the input and state predictions for the entire forecast horizon, N. So the input u ∈ R6N in J1 (u) and J3 (u) contains the N stacked input vectors u ∈ R6 for each time step in the prediction horizon and is expressed as T u = u(k |k )T , u(k + 1|k )T , ..., u(k + N |k )T . Similarly, x ∈ R6N contains the predicted T state variables for the entire horizon, x = x (k |k )T , x (k + 1|k )T , ..., x (k + N |k )T which 35 are inherently functions of u. The forecasted state variables are calculated using Equation 3.9 using the predicted input variables, u, the temperatures at the current time step, x̂ (k ) ∈ R6 , and the known stacked outdoor temperature over the entire horizon, T w(k |k )T , w(k + 1|k )T , ..., w(k + N |k )T . The setpoint temperatures are considered as constant values and are the same for every zone except for the attic, which isn’t thermally regulated. In J2 (u), Tsp ∈ R6N is the column vector [22, 22, 22, 22, 22, ∞] T stacked N times for the case of a constant temperature setpoint. When the variable temperature band is used, Tsp is taken as some offset from the heating setpoint at the current time step, represented as Tsp = [ Tmin , Tmin , Tmin , Tmin , Tmin , ∞] T + σ. Where σ is used to direct the temperature tracking goal to be closer to the center of the temperature band as opposed to following the heating setpoint. The constraints in Equation 3.8 are applied only to the five zones that are thermally regulated and are set to the values shown below. u = 5000 W, x = Tmin , x = Tmax (3.16) where Tmin and Tmax are the heating and cooling setpoints, respectively, that are defined in Table 3.4. This optimization problem formulation and MPC control algorithm were implemented in MATLAB R2018a with a basic academic license. To solve the optimization problem, CVX, a package for specifying and solving convex programs was used [67]. The optimization problem was formulated as a linearly constrained quadratic optimization problem in order to satisfy the disciplined convex programming rules that constrict the solution feasibility in CVX. 3.4 Results Three separate studies were performed in order to present a cross-comparison between different control methodologies and between the same MPC control scheme utilizing different parameters. The first is a comparison between the basic on/off control scheme and the MPC scheme implemented with on/off heating control and variable heating control. Next, a brief look into the impact of small variations in the forecasting horizon is presented. Finally, a parametric analysis based on the weighting of the third cost function, J3 (u), 36 is shown by varying φ across a range that depicts the effect of this control parameter in eliminating demand spikes and leveling temperature fluctuations. 3.4.1 Comparison to Basic On/Off Control In this section, results for the simulations performed using the three different control techniques are shown for comparison. They are compared numerically based on the total energy usage over a two-day period, which is calculated by integrating the power input over the total 48-hour simulation period. The total heating energy usage over a 48-hour period in January is shown in Table 3.8. In this case, the lowest energy consumption is obtained using the basic on/off controller. This is because the thermal comfort cost was weighted more heavily than the heat input cost for the MPC control, which limits its usefulness as an energy reduction strategy because the main priority of the controller is the setpoint tracking. A time step of 60 seconds was used for the forecasting in this section, which makes the optimization problem have a large computational cost for large forecasting horizons. In order to maintain a subhour simulation time, only forecasting horizons of less than 30 minutes are considered in this section. The forecasting horizon of MPC control in building thermal systems is an important parameter because it limits the disturbance horizon that the controller considers when setting control laws in real time. This future disturbance prediction is one of the main benefits to an MPC controller, as it allows the system to utilize passive components to overheat or overcool the building to minimize the total energy consumption over some time frame [68]. Seeing that a horizon of less than 30 minutes doesn’t capture significant temperature fluctuations, the energy usage of the MPC control laws is not expected to cause significant reductions in energy usage. Accordingly, the implementation of MPC is presented here as a reference case study that provides an example of how MPC can be implemented with a simple building model. Table 3.8. Total heating energy usage for 48-hour period in January five different control simulations. MPC10, MPC20 and MPC30 are MPC variable control simulations using a forecasting horizon of 10, 20 and 30, respectively. Basic On/Off 1132.44 Total Heating Energy (MJ) MPC On/Off MPC10 MPC20 1150.20 1157.05 1152.24 MPC30 1151.25 37 The temperature response and heat input are shown for the three different control techniques in Figure 3.4, Figure 3.5 and Figure 3.6 for a simulation period of 48 hours and one hour. For this comparison, a forecast horizon of 20 minutes is used. There are three major distinctions in the behavior of these three control techniques. First, both of the on/off control methods show unstable behavior that could be problematic in cases where high-frequency demand changes are a concern either for the internal building systems or the distribution network. These figures show that on/off cycles are occurring between sequential time steps of 60 seconds. This is often an undesirable situation for real building system operation and indicates that the heat input for this simplified thermal model of the small office building could be oversized. The second distinction is in the difference in setpoint tracking trajectories for each of the five zones. In the basic on/off control, the zone temperature control is isolated for each of the separate zones, while in MPC control all of the zones are considered simultaneously in the computation of the optimal control input. The third distinction is in the ramp-up behavior, which is shown in the single hour graphs in Figure 3.4, Figure 3.5 and Figure 3.6. The basic on/off controller immediately turns on the heat input to all zones for the same reason as what is seen in distinction 2, while the MPC controllers first activate the heat supply for the external zones. This is a direct result of the model-based control technique. In order to test the hypothesis that decreasing the forecasting horizon decreases the ability of the MPC controller to minimize energy usage, a parametric study was conducted by running a 48-hour simulation with a horizon of 10 minutes, 20 minutes and 30 minutes. The results for this are shown in Figure 3.7 and the difference in total energy consumption for the three different forecast horizons can be seen in Table 3.8. We see that the total heating energy is reduced as we increase the forecasting horizon, which is an effect of providing a larger domain for CVX to obtain an optimal solution. The difference in energy consumption by changing the horizon from 10 to 30 minutes is only 0.5%. However, we expect that this change would be much more favorable if we increased the horizon to exceed one hour because this allows the MPC controller to account for larger temperature fluctuations when evaluating the objective functions. Results for a similar MPC simulation using a more ideal forecasting horizon are shown in Section 3.4.2. 38 (a) (b) Figure 3.4. Control input and state response trajectory for basic on/off control over a 48-hour period (a) and for a one-hour period (b) in January. 39 (a) (b) Figure 3.5. Control input and state response trajectory for MPC on/off control over a 48-hour period (a) and for a one-hour period (b) in January. 40 (a) (b) Figure 3.6. Control input and state response trajectory for MPC variable control over a 48-hour period (a) and for a one-hour period (b) in January. 41 Figure 3.7. Forecast horizon comparison for MPC variable control simulation with a forecast horizon of 10, 20 and 30 minutes. 3.4.2 Demand Peak Shaving One of the highly noted benefits of MPC is its ability to support demand response strategies. Load patterns over a day are typically easy to predict for an office scenario, so the application here is a relatively trivial case. However, when buildings are interconnected with the electrical distribution system or in some cases with other buildings in the case of a campus, there are a multitude of opportunities to shift or reduce demand in a beneficial manner. Benefits associated with increased load flexibility include cost savings, potential emissions or water usage reductions, aiding grid operations and reliability, and providing opportunities for total energy savings. The results provided here demonstrate a simple case where occupancy forecasts can be used to reduce the ramping burden and the associated demand spike that comes with switching between thermal comfort setpoints. For the simulations performed in this section, a time step of 15 minutes is used to reduce 42 the computational expense of using a 24-hour prediction horizon. This time step proved to be sufficient in relation to the thermal response of the building based on a comparison with the results that were obtained using a one-minute time step. The temperature trajectory for these MPC simulations using three different values of φ is shown in Figure 3.8. In each of these graphs, we see that the main effect of changing the weight on the heat input time derivative objective occurs in the periods of changing temperature setpoints, which is better explained in Figure 3.9 where the heat input throughout the simulation is plotted. The temperature profile shows that the cost associated with deviation from the heating setpoint offset used in J2 (u) is outweighed by the heat input cost, since the temperature follows the lower constraint, Tmin , instead of the offset temperature setpoint, Tmin + σ. While this is sufficient in terms of thermal comfort, it would be more desirable to have the temperature track a setpoint more towards the center of the band as it does for the constant setpoint simulations. This could be accomplished by changing the weights of the objective functions, J1 (u) and J2 (u). Figure 3.8. Comparison of temperature trajectories for MPC control utilizing a difference cost weighted with φ = 0, φ = 1000 and φ = 10000. 43 Figure 3.9. Comparison of heat supply trajectories for MPC control utilizing a difference cost weighted with φ = 0, φ = 1000 and φ = 10000. 3.5 Summary The simulations presented here give insight into how an MPC temperature tracking controller operates in relation to basic on/off control and how its operation can be adjusted by tuning the weights in the objective function. All of the control simulations operated on a simple RC thermal circuit model that was constructed based on the geometry of a reference small office building. This crude modeling approach proved to be useful in mimicking the thermal response that would be expected from a real building and therefore has useful applications in developing and testing building energy control algorithms, as is shown here and in other instances [61], [62], [63], [64]. The three different control techniques shown in Section 3.4.1 were designed in a way that successfully maintained the temperature sufficiently close to the setpoint. However, the unstable behavior of the heat input for the on/off control methods indicates that these algorithms should be modified to reduce the occurrence of high-frequency changes that could impact thermal comfort, building equipment performance, and power demand requirements. This could potentially be accomplished either by reducing the constant heat 44 input (i.e., resizing the heating system) or adding control logic to the algorithm for on/off control. The opportunity for demand management with the utilization of MPC is extensive, as it allows adjustments for load shifting and peak shaving to be automated. One technique for peak reduction was shown here by simply adjusting the weight distribution in the objective function. A similar method could be used by a building energy manager in real time as priorities and operative requirements evolve. The advantage of this flexibility also extends to situations where energy demand is met by microgrids and district heating infrastructures. In these cases, which are increasingly utilizing intermittent renewable electricity generation, demand flexibility can significantly reduce the supply side burden. CHAPTER 4 UNCERTAINTY IN BUILDING PERFORMANCE SIMULATION The process of modeling a building, and the systems that sustain its operation, has been an essential part of the design and maintenance of buildings for the majority of the 20th century. In particular, building performance simulation (BPS), the process of using computational mathematical modeling to represent building performance characteristics, has been in a state of rapid development and improvement since the late 1960s [69]. Methods have progressed from crude models that neglect many aspects of the internal physics of building energy systems, such as the one presented in Chapter 3, to advanced simulations that account for finely detailed processes in the modeled building. This progression through the 1990s is described succinctly by Clarke and Maver and progression through the present day is summarized in Table 4.1 [69], [70]. There is a wide array of software packages that are available to assist with the process of BPS; 20 of the most well-known software packages are described and compared by Crawley et al. [11]. In addition to improvements in the physical representation of buildings, the application of BPS has extended beyond the initial process of building design and commissioning and is now acting as an essential managing agent throughout the entire life cycle of buildings. The expansion of the domain in which BPS is used has introduced new challenges and highlighted existing challenges that have been the focus of improvement in recent years (Table 4.1). Hong et al. have listed 10 challenges that they expect to be central to new research and development in the field of building simulation [13]. Many of these challenges have been in the scope of building simulation research in the past decade but remain to be major technical challenges, such as addressing the gap between predicted and actual building performance and modeling human occupant behavior. These things not only 46 Table 4.1. Building performance simulation advancement timeline for the 1960s through present day. 1st Generation (1960s and early 1970s) [70] Handbook oriented, simplified, piecemeal. 2nd Generation (1970s) [70] Dynamics important, less simplified, still piecemeal. 3rd Generation (1980s) [70] Field problem approach, numerical methods, integrated energy subsystems, heat and mass transfer considered, better user interface. 4th Generation (1990s) [70] Computer aided design integration, advanced numerical methods, intelligent knowledge base, advanced software engineering. 5th Generation (2000s) User friendly and open source (EnergyPlus), inclusion of occupancy behavior and comfort, spanning multiple technical domains. 6th Generation (2010s) Improved visualization capabilities, building work flow and life cycle integration, advancing capabilities for application in real time operation and control. Indicative, applications limited, difficult to use, computers not available for BPS applications. Broad ranging applications, available for most skill levels, open source, spans multiple technical domains, fully integrated with computers and computer science. 47 affect the model’s accuracy at face value, but they also impact the usefulness of BPS in continuous building operation and controls. For BPS to be used with confidence in live building operation and control, the model must closely represent the actual behavior of the building [71]. But the resources, skill-level, and time required to calibrate building models to the extent that is necessary are expensive and in some cases impractical. As a consequence, BPS has not been adopted on a large scale in building control systems [58]. A typical building model calibration process consists of perturbing select model input parameters in order to minimize the residual. It has been noted that this process can be problematic because selecting the correct input parameters to vary and the magnitude in which they are varied strongly depends on the expertise of the building modeler [72]. This is especially an issue today, since building models are becoming increasingly complex and often contain thousands of input parameters. Performing an uncertainty analysis to relate uncertainty in different model inputs to the corresponding impacts on uncertainty in the model output can provide insight into the inner workings of the model and assist in parameter space reduction [70]. And, as a more refined version of uncertainty analysis, a sensitivity analysis can directly point the designer to the most influential inputs [73], [74]. Approached in this sense as a parameter reduction problem, uncertainty and sensitivity analysis done on the building model can significantly reduce the calibration burden by singling out the model parameters that impact the desired outputs. 4.1 Problem Definition This work intends to investigate the usefulness of a procedure that can point a building modeler to parameters that should be the focus of calibration efforts and to identify parameters that should be considered for control variables in system control schemes. A procedure that has been established in the BPS research community [70], [75] is applied to a simple building model to identify the input parameters that influence output uncertainty for electric energy and thermal comfort calculations. For simplification purposes, EnergyPlus input parameters are first grouped into three different categories: schedules/setpoints, materials, and equipment. This allows the uncertainty analysis results to be viewed categorically and acts as a high-level filter for identifying high-impact parameters. A global sensitivity analysis is then performed to extract the individual parameters that have 48 the highest impact on output uncertainty. This is done separately for electric energy and thermal comfort outputs. The results for this sensitivity analysis are compared directly to the results from the uncertainty analysis in order to confirm that the parameters identified through the uncertainty analysis are, as expected, responsible for most of the uncertainty in each of the three categories. This is performed as a precursor for a study into the application of EnergyPlus building modeling for use in MPC for building energy system controls. In order to gain insight into how uncertainty in high impact parameters affects the performance of an MPC scheme that utilizes EnergyPlus building models, we briefly look into how uncertainty impacts the building performance optimization portion of MPC. This is done by running an optimization scheme in an MPC utilization framework with uncertainty applied to parameters that are identified through the sensitivity analysis. 4.1.1 System Model The EnergyPlus model that is being used in this study is a small office reference building from the DOE commercial prototype building model database [60]. This is the same building that was modeled in Chapter 3. The building envelope and interior are shown in Figure 4.1. The building has five conditioned thermal zones that interact with each other through internal walls, four of which are external and exchange heat with the surrounding environment through external walls. In addition to heat exchange with spaces horizontally, each thermal zone exchanges heat through the ground and with the attic of the building. Figure 4.1. Isometric (left) and top view of the interior (right) of the small office prototype building used in this study. 49 In this study, building performance is being quantified using two different metrics: facility electric energy (FEE) and predicted mean vote (PMV) [76]. PMV is a thermal comfort metric and here it is calculated using the Fanger comfort model, which is the oldest and most established method for quantifying thermal comfort [77], [78]. The Fanger thermal comfort model is constructed based on the modes of heat transfer that affect the heat exchange between a person and their environment. The detailed formulation of the method is included in [78]. PMV is a discrete index that ranges from -4 to 4, describing the different levels of sensation listed in Table 4.2. While PMV is defined on a discrete scale, EnergyPlus reports the calculations on a continuous scale. A building simulation period of two weeks is used in order to sufficiently capture the behavior of the building through weather and occupancy changes. Using a simulation period of only two weeks also significantly reduces simulation run time and speeds up the uncertainty and sensitivity analysis. The location that is used for the results presented here is Lansing, MI, which is in IECC climate zone 5A. 4.1.1.1 Nominal Performance The nominal performance outputs are defined as the outputs that result from using the predefined input parameter values that are contained in the small office prototype model. The simulations here are initially studied over a summer two-week period (July 16 - July 29) and a winter two-week period (January 15 - January 28), but some of the results are only presented for the summer period because seasonal differences aren’t the focus of this study. The nominal outputs for the summed FEE and averaged PMV values over these two-week intervals are given in Table 4.3. While there is substantial variance between the PMV values in different thermal zones, this study uses the zone average calculations for evaluating thermal comfort. The hourly values for the building performance over the summer and winter simulation periods are shown in Figure 4.2. Table 4.2. Predicted mean vote sensation scale used for thermal comfort quantification. very cold -4 cold cool -3 -2 slightly neutral slightly warm cool warm -1 0 1 2 hot 3 very hot 4 50 Table 4.3. Nominal building performance metrics calculated over two-weeks in the summer and winter. Two-week sum (facility electric energy) and two-week average (predicted mean vote). Summer Winter Zone 1 PMV -0.094 -1.053 Zone 2 PMV -0.067 -1.169 Zone 3 PMV -0.148 -1.168 Zone 4 PMV 0.019 -1.169 Core Zone PMV -0.043 -0.672 Zone Ave PMV -0.066 -1.046 FEE (MJ) 6792 6845 Figure 4.2. Nominal performance for small office building simulations over summer (07/16-07/30) and winter (01/15-01/29) simulation periods. 4.1.1.2 Input Parameters The simulated uncertainty in this study is initiated by varying the values of the EnergyPlus input data file numerical parameters. A total of 907 parameters are included in the analysis and are broken down into three categorical groups: schedules/setpoints (221 parameters), materials (576 parameters), and equipment (110 parameters). The input data file class types that were included in each of these groups are listed in Table 4.4. All of the numerical parameters for each class were included in the uncertainty and sensitivity analysis. 51 Table 4.4. EnergyPlus IDF classes included in UA/SA. Schedules/Setpoints (221 parameters) SCHEDULE:COMPACT, SETPOINTMANAGER:SINGLEZONE:REHEAT, PEOPLE Materials (576 parameters) WINDOWMATERIAL:GLAZING, MATERIAL, MATERIAL:NOMASS, WINDOWMATERIAL:GAS Equipment (110 parameters) LIGHTS, PLANTLOOP, AVAILABILITYMANAGER:NIGHTCYCLE, WATERUSE:EQUIPMENT, ELECTRICEQUIPMENT, EXTERIOR:LIGHTS, DESIGNSPECIFICATION:OUTDOORAIR, WATERHEATER:MIXED, FAN:ONOFF, COIL:COOLING:DX:SINGLESPEED, COIL:HEATING:DX:SINGLESPEED, COIL:HEATING:GAS, CONTROLLER:OUTDOORAIR, AIRLOOPHVAC:UNITARYHEATPUMP:AIRTOAIR, PUMP:CONSTANTSPEED 4.2 Uncertainty and Sensitivity Calculation The three different categories (schedules/setpoints, materials, and equipment) were selected and constructed based on personal judgment to make the analysis intuitive with respect to the physical system hierarchy (i.e., the classes included in each category listed in Table 4.4 were chosen by the researcher). Here, the uncertainty analysis (UA) is presented first as a preliminary investigation into the weight that the three different categories of parameters have on the model output uncertainty. The UA was done by applying uncertainty distributions to each of the input parameters individually. For the majority of the parameters, a 20% uniformly distributed uncertainty interval was applied to the nominal value. However, there were cases where this was either too extreme for the simulation constraints or wasn’t effective. For example, the heating and cooling temperature set points could not be varied by 20% because it would cause the heating setpoint to overlap with the cooling setpoint, a case that is highly problematic in HVAC control. Another example of this was any object containing radiation properties (absorptive, reflectively and/or transitivity) 52 because these parameters have interdependent constraints forcing their sum to be less than unity. Accounting for these anomalies, the logic that was used for assigning uncertainty is described by the pseudocode in Table 4.5. Additional rules and reasoning could be applied to make these uncertainty intervals more reasonable within the contextual building system. But due to the large amount of parameters and limited time, a homogeneous uncertainty interval was applied with the few exceptions that are mentioned here. In order to make the problem of calculating sensitivity indices for the large parameter set that is used here more tractable, a Gaussian Process Regression (GPR) model is fit to the EnergyPlus model for the simulation periods and climate regions specified. GPR is a probabilistic method and works well with multivariate sets that follow a Gaussian distribution as is the case with building energy models [79]. GPR is generally robust to hyperparameters and fitting kernels, which reduces the time required for optimizing the regressor to individual problem characteristics [80]. In the context of building energy modeling, GPR has been shown to perform similarly to or better than other regression techniques such as support vector regression, kernel ridge regression, and neural networks Table 4.5. Brief of algorithm used to apply simulated uncertainty intervals to input parameters. Result: Returns a unique uncertainty interval for all parameters in defined set. for each parameter in set do uncertainty interval = 20%; extract nominal value; extract acceptable range; if temperature setpoint parameter then reduce uncertainty interval to 5%; else if solar radiation, lighting, or glazing parameter then constrain corresponding α, ρ, τ to be less than 1; else if nominal value is zero then use right triangular distribution between zero and 10% of the nominal value; else apply 20% uniform uncertainty distribution; end end 53 [80], [81], [82]. However, this is highly dependent on the problem type. The fitting of a metamodel allows us to avoid iterative EnergyPlus simulations that have a high computational cost and also provides a continuous function that supports the use of direct derivative based methods for calculating global sensitivity indices. One downside to fitting a metamodel to an EnergyPlus simulation is that it reduces the applicable model domain. Because the metamodel is fit to the EnergyPlus simulation outputs over a specific time frame and with a particular set of weather data, the accuracy of the model typically diminishes outside of the domain spanned by the training data (i.e., the time period and location). This makes the metamodeling approach only useful when a live or responsive output isn’t a requirement. The metamodeling technique described here is commonly used for building optimization, parameter reduction and calibration [83], [74], [75], [84], [85]. A GPR technique was selected because, in addition to reasons mentioned previously, it has been shown to perform well when fitting the high dimensional and nonlinear data sets that are posed by building energy systems [84], [86]. The GP regressors were applied in Python using the scikit-learn Gaussian processes toolbox [87], [88]. This toolbox provides an extensive platform supporting many different machine learning techniques. It has a host of kernels that are available for fitting training data and has functionality that makes it easy for researchers to search for optimal hyperparameters and kernels. Here, a superposition of DotProduct and WhiteNoise kernels was used to fit the model to EnergyPlus data for the specified simulation spatial and temporal domain. Once the metamodels were fit using GPR for both the FEE and PMV outputs, a sensitivity analysis (SA) was done on all 907 parameters. The SA analysis is expected to support the UA analysis results since the majority of the high impact parameters found in the SA analysis should belong to the high impact categories that are found in the UA analysis. SA has been used extensively within the BPS research community for screening, model calibration, and optimization [75], [89], [90]. As indicated by Eisenhower et al. and Zheng et al., variance-based sensitivity methods do not always capture the full extent of behavior in energy models, making derivative based sensitivities a potentially better option [83], [75]. Based on this assertion, the Derivative-based Global Sensitivity Measure (DGSM) presented by Sobol and Kuchereneko [91] is used here . The sensitivity indices are evaluated using 54 the SALib DGSM module for Python [92]. The DGSM sensitivity index corresponding to the facility total electric energy usage function, f 1 ( x0 , x1 , ..., x906 ) and the zone average PMV function, f 2 ( x0 , x1 , ..., x906 ) for each parameter, xi for i = {0, 1, 2, ..., 906}, is tot Si,j = ( ximax − ximin )2 π2 Dj Z Hn ∂ fj ∂xi 2 dx (4.1) where ximax and ximin are the max and min acceptable values, respectively, for parameter i, D j is the total variance of the output function f j ( x0 , x1 , ..., x906 ) and H n is the sample space formed by the intersection of the bounds on the input parameters, where n = 907. The DGSM indices are evaluated using the SALib DGSM module by numerically computing the integral in Equation 4.1 over a large set of randomly sampled points in H n . The usage of the DGSM index assumes that the output function is differentiable over the entire domain, H n [91]. This is satisfied here, since the functions, f 1 ( x0 , x1 , ..., x906 ) and f 2 ( x0 , x1 , ..., x906 ), that are fit using GPR are infinitely differentiable [93]. It has also been noted that for highly nonlinear models, the DGSM index ranking procedure can produce false conclusions [91], [94]. A GPR model, f (x) = h(x)T y, is a linear combination of the training target values, y, weighted with some function, h(x), of the inputs, x, (i.e., a linear smoother) [93]. This is different from a linear model where the output is a linear combination of the inputs, x. The DGSM index is still used here to take advantage of its computational efficiency and based on the similar use of derivative-based sensitivity indices in past studies [74], [83]. However, an additional level of scrutiny is required when ranking parameters with the DGSM index because of its unsurety for nonlinear models. This is done here by justifying high-ranking parameters based on their analytical or pseudo-physical relevance in the model. 4.3 Building Performance Optimization Building performance optimization is a major functional aspect of MPC. In order to better understand the effect of uncertainty in this layer of MPC algorithms, a simple building performance optimization procedure is presented. The optimization problem formulated here is applied for five different uncertainty scenarios that come directly from the sensitivity analysis previously explained. The optimization problem formulated here is for two variables, a constant heating setpoint, Tmin , and a constant cooling setpoint, Tmax . By using constant heating and cooling 55 setpoints, the complexity of the optimization problem is significantly reduced from the ideal case where the optimization algorithm would optimize the setpoints on an hourly basis. However, in this case, it is sufficient and perhaps more effective for showing the effects of uncertainty on optimized performance predictions. A multiobjective optimization approach was used here, so instead of assigning weights to different objective functions and obtaining a single optimal solution, all feasible solutions that minimize each objective individually is found and presented as the Pareto frontier. A Pareto frontier represents the set of points that cannot be improved further in one objective without degrading the other objective. The problem is formulated using the two decision variables, Tmin and Tmax , two objective functions, J1 ( Tmin , Tmax ) and J2 ( Tmin , Tmax ), and 48 constraints. min (J1 ( Tmin , Tmax ), J2 ( Tmin , Tmax )) s.t. −1 ≤ f pmv (t) ≤ 1 for t = 1, 2, ..., 48 (4.2) with the two objectives being J1 ( Tmin , Tmax ) = f et Q f e (4.3) t J2 ( Tmin , Tmax ) = f pmv R f pmv (4.4) where f e ∈ R48 and f pmv ∈ R48 are the EnergyPlus HVAC energy usage and PMV calculations at hourly intervals over the entire 48-hour simulation period, respectively. Here, Q ∈ R48×48 , R ∈ R48×48 are scaling matrices, as opposed to weighting matrices as used in Chapter 3 that are used to evenly scale the two objectives. They are defined and utilized in a similar manner as they are in Chapter 3 but fundamentally have a different effect since this is solved as a multiobjective optimization problem. The weighting matrices are calculated as Q = ρI, R = ψI (4.5) where ρ = 1e-7 and ψ = 5 were found by adjusting their values until both of the objectives were similar in magnitude. This optimization problem was implemented using Platypus, a Python framework for evolutionary computing that specializes in multiobjective evolutionary algorithms [95]. Within the Platypus framework, a very popular multiobjective genetic algorithm, NSGA-II, was used. Those wishing to know more about NSGA-II can reference [96]. One hundred problem evaluations were used for performing every optimization presented here. In order 56 to construct the optimization framework in a way that allows the optimization algorithm to treat an EnergyPlus building model as an objective function, a Python EnergyPlus co-simulation platform, eppy, was used [97]. This made the process of updating EnergyPlus input data files and iterating through EnergyPlus simulations relatively straightforward. 4.4 Results The results are first presented for the uncertainty analysis (UA), showing a crude breakdown for the identification of high impact input parameter types. The effect of input uncertainty in the three categories is quantified for the two simulation outputs, FEE and zone averaged PMV. Other metrics that are commonly considered in the design and operation of buildings are also calculated as a reference for those interested. Namely, primary energy consumption, electricity cost, and indirect emissions, which are all assumed to be proportional to electric energy usage for this analysis. The results for the sensitivity analysis (SA) are presented as a follow up to the uncertainty analysis. The SA results are compared to the UA results as a point of discussion in terms of validation and verification. They are also used directly in Section 4.4.4 when we present the effects of having uncertainty in the most sensitive parameters on the optimization results in the simplified building performance optimization problem being used here. The output distributions are presented for the UA and SA in the form of box plots. In these, data are split up by quartiles and represented as is shown in Figure 4.3. Figure 4.3. Interpretation of box plots that are used to show output distributions. 57 4.4.1 Uncertainty Analysis By varying the inputs to the EnergyPlus model by 20%, the results are not expected to be trivial because the model complexity and interactions between all three of the categories that are described here amplify noise in the complex building system. Splitting the input parameters into the three different categories aforementioned allows us to increase our perspective on how these different model components interact with each other and ultimately affect the model output. Here, the model was run with 100 samples for each category, and the standard deviation in the output distribution was used as a crude measurement for uncertainty. Figure 4.4 shows a daily profile for FEE and PMV with the 95% confidence interval for each hour calculated with 100 simulation samples. Figure 4.5 shows the output uncertainty distributions for aggregated output metrics; the two-week summation for FEE and the two-week average for PMV. Both of these figures, although more easily interpreted in Figure 4.5, show that for FEE, schedules/set points and equipment parameters are more influential in the output uncertainty compared with the 576 material parameters. A similar understanding can be drawn from Figure 4.5 for PMV; however, for this output metric, the equipment parameters play a minor role. Figure 4.4. Daily trends and 95% confidence interval relating to each category that 20% simulated uncertainty was applied to separately. 58 Figure 4.5. Output distributions for simulated 20% uncertainty the three separate categories and for results related to simulating uncertainty in all parameters. Total facility electric energy is an output that can be used to compute values that are commonly of interest to building designers, operators, and occupants. Here, relevant factors are used to convert from facility electricity usage to emissions, coal usage, and electricity cost for the two-week simulation period. The results for these calculations are shown in Table 4.6. These are meant to provide insight into how uncertainty in the input can not only impact the direct model output but also many other metrics that are influential in decision making processes throughout the lifetime of the building. The factors that are used for the emissions calculations and the coal usage estimates were computed using a software tool that provides estimates for the fuel mix and associated emissions at hourly time intervals for regions within the United States as described in Chapter 2 [98]. A transmission/distribution loss factor of 6.5% and a fuel conversion efficiency factor of 30% was applied for the computation of total coal usage. Note that the estimate for coal usage uses a fraction for the amount of electricity that is being supplied by coal out of the entire generation pool for each hour, so this is not assuming that all of the electricity is supplied by coal. A unit price for electricity in Lansing, MI was taken to be 10.93 cents/kWh from Electricity Local [99]. 59 Table 4.6. Uncertainty analysis for primary energy, emissions and cost quantification over a two-week summer simulation period. mean, µ std, σ CO2 [kg] 1192.6 45.829 NOx [kg] 0.8275 0.0318 4.4.2 SO2 [kg] 1.2619 0.0485 coal consumption [kg] 54873 38823 cost [$] 205.09 7.8928 Sensitivity Analysis The results for this sensitivity analysis are presented as a way to narrow down the parameters within the scope of a building model to simplify the perceived input-output relationships. To reduce the computational burden of doing a global sensitivity analysis for a large input parameter set, a Gaussian Process Regression (GPR) model was trained to reflect the behavior of the EnergyPlus model for a specified temporal and spatial domain. This provides a closed functional form that can be used in the computation of DGSM indices, as discussed in Section 4.2. A convergence analysis was performed by training models with parameter sets ranging from 100 to 5000 samples, and the Bhattacharyya distance was calculated to identify the point at which the GPR reached a stable solution. The Bhattacharyya distance is a metric that quantifies the statistical variation between two normal distributions [100]. Referencing Figure 4.6, the statistical variation between the EnergyPlus and GPR models converges at a sample size of approximately 4000 samples. Hence, from this point forward the GPR that was trained with 4000 samples is utilized. Figure 4.7 presents a comparison of the GPR model output distributions to the EnergyPlus output distributions with 5000 samples and Table 4.7 provides a statistical comparison for the two models. While the PMV GPR fit showed a close statistical match, the FEE GPR model wasn’t able to match the standard deviation of the EnergyPlus model. We expect that this is a result of applying an uncertainty distribution to a highly sensitive parameter that shouldn’t necessarily be considered in the analysis. This will be more closely analyzed in the future so that the source of error can be identified. For now, the fits are considered sufficiently accurate for the purpose of identifying a set of model input parameters that have a high DGSM index. The DGSM indices that were calculated using the GPR model are shown in Figure 4.8 and are categorized into the three groups that are used in this chapter for comparison. The 60 Figure 4.6. Gaussian process model convergence analysis for two-week summer simulation. Bhattacharyya distance as a function of training sample size is plotted to represent the convergence of EnergyPlus (EP) and Gaussian Process (GP) model output distributions for 5000 samples. Figure 4.7. Comparison of EnergyPlus model with Gaussian process model that was trained with 4000 simulation samples. Output distributions are for a sample size of 5000 points. 61 Table 4.7. Statistical comparison of output distributions for EnergyPlus model and Gaussian process model for 4000 training samples. EnergyPlus GPR error % µ electric [MJ] 6741.0 6741.5 0.0073 µ PMV -0.2085 -0.2088 0.1240 σ electric [MJ] 265.58 110.10 -58.544 σ PMV 0.4002 0.3849 -3.8163 Figure 4.8. Derivative-based Global Sensitivity indices for facility electric energy output and predicted mean vote output. All indices are included in the left two subplots and the top 10 parameters for both output variables are included in the right two subplots. 62 10 input parameters with the highest DGSM indices were selected for both model outputs and are shown in the two right plots of Figure 4.8. These 10 input parameters are ranked with their category, name, and DGSM index value in Table 4.8 and Table 4.9 for FEE and PMV, respectively. The sensitivity analysis results are presented only for summer simulation periods because winter simulations yielded nearly identical results. 4.4.3 Uncertainty and Sensitivity Analysis Comparison In Section 4.4.2, we found that the sensitivity of the top 10 parameters for FEE and PMV outputs were significantly greater than the sensitivity of the remaining parameters. In this case, it is expected that a reduced order GPR model trained only using these high impact parameters will maintain the same relation to the original EnergyPlus model as the original GPR model that was trained using all 907 input parameters. This is verified in Figure 4.9 where the original EnergyPlus and GPR models are compared to reduced-order GPR models that were trained using only the top 10 and top 5 parameters for both outputs. The statistical comparison of these reduced-order models to the EnergyPlus model with uncertainty applied to all input parameters is given in Table 4.10. We see that the FEE standard deviation for both reduced models has a large error of close to 60%, which is close to the 58.5% error for the full order GPR model. In addition to this, an interesting result is shown in Table 4.11 where the reduced order GPR models are compared Table 4.8. Top 10 parameters, based on DGSM index, for facility electric energy output. Rank ID 1 2 3 4 5 6 7 8 9 853 839 880 865 870 860 875 653 188 10 730 1 Off Category Equipment Equipment Equipment Equipment Equipment Equipment Equipment Materials Schedules/ Setpoints Equipment Parameter Name Heater1 SHWSys1 Water Exterior Lights B PSZ-AC:5 Fan PSZ-AC:2 Fan PSZ-AC:3 Fan PSZ-AC:1 Fan PSZ-AC:4 Fan Std Wood 6inch Activity Schedule DGSM Index 0.9505 0.1809 0.0801 0.0772 0.0431 0.0344 0.0291 0.0203 0.0198 200mm Normalweight con- 0.0158 crete floor Cycle Parasitic Fuel Consumption Rate 63 Table 4.9. Top 10 parameters, based on DGSM index, for predicted mean vote output. Rank 1 ID 188 2 192 3 189 4 190 5 6 7 740 343 18 8 9 10 853 741 854 Category Schedules/ Setpoints Schedules/ Setpoints Schedules/ Setpoints Schedules/ Setpoints Materials Materials Schedules/ Setpoints Equipment Materials Equipment Parameter Name Activity Schedule DGSM Index 0.9527 Clothing Schedule 0.1147 Work Eff Schedule 0.0109 Air Vel Schedule 0.0075 CP02 CARPET PAD (Rth ) LoE TINT 6MM Bldg Equip Schedule 0.0018 0.0011 0.0005 SHWSys1 Water Heater1 CP02 CARPET PAD (α) SHWSys1 Water Heater2 0.0004 0.0004 0.0003 1 Off Cycle Parasitic Fuel Consumption Rate 2 Off Cycle Parasitic Heat Fraction to Tank Figure 4.9. Comparison of EnergyPlus model with Gaussian process models that were trained with all 907 parameters, top 10 parameters and top 5 parameters based on the DGSM index ranking. 64 Table 4.10. Statistical comparison between EnergyPlus simulations with 20% uncertainty applied to all parameters and EnergyPlus simulations with 20% uncertainty applied to the top 10 parameters and top 5 parameters for FEE and PMV using a sample size of 5000. Top 10 Baseline Reduced error % Top 5 Baseline Reduced error % µ FEE [MJ] µ PMV σ FEE [MJ] σ PMV 6741.0 6791.9 7.5508 -0.2085 -0.2331 11.799 265.58 102.91 -61.159 0.4002 0.3997 -0.1249 6741.0 6792.4 7.6249 -0.2085 -0.2365 13.429 265.58 98.867 -62.773 0.4002 0.4075 1.8240 Table 4.11. Statistical comparison of output distributions for EnergyPlus model and Gaussian Process Regression (GPR) model for 4000 training samples using only the top 10 and top 5 parameters identified using DGSM sensitivity indices (4.8, 4.9). Baseline EnergyPlus Top 10 EnergyPlus GPR error % Top 5 EnergyPlus GPR error % µ FEE [MJ] µ PMV σ FEE [MJ] σ PMV 6741.0 -0.2085 265.58 0.4002 6791.9 6791.9 4.73e-07 -0.2331 -0.2323 -0.0030 102.91 102.90 -4.19e-05 0.3997 0.3740 -0.0644 6792.4 6792.4 -1.83e-07 -0.2365 -0.2355 -0.0044 98.867 98.8648 -2.49e-05 0.4075 0.3806 -0.0661 to EnergyPlus simulations with uncertainty only in the top 10 and top 5 parameters. The fact that the standard deviation error is almost completely eliminated when we only apply uncertainty to the high impact parameters when running the EnergyPlus simulation further supports the hypothesis that the EnergyPlus simulations are being skewed by unstable input parameters, as suggested earlier. It is also possible that important parameters were not weighted properly in the training process of the full order GPR models and, in turn, were eliminated from the reduced order models. This is an unfortunate consequence of poor fitting in the field of machine learning. The results obtained here are still impactful, however, since the parameters identified here are still considered highly sensitive in comparison to the rest of the input set. 65 Figure 4.10 is a final comparison for the aggregated output calculations that were shown previously in Figure 4.5 to the same metrics calculated from EnergyPlus simulations that had uncertainty applied only to the top 10 and top 5 parameters. We see that for the PMV output, the majority of the uncertainty is still captured in the top 10 and top 5 distributions. Alternatively, it is apparent that some of the equipment parameters that account for a large portion of the uncertainty in the FEE output were not included in the top 10 and top 5 parameter sets. This highlights the same GPR fitting issue that was discussed previously. 4.4.4 Uncertainty Impact on Optimization In an MPC framework, building performance optimization is ingrained into an automated procedure where the optimal control inputs are calculated and implemented recurrently on various time scales, for example, once every minute or 15 minutes as was seen in Chapter 3 or once every hour as would be the case here if the optimization problem was adopted in an MPC scheme. For this reason, it is a vital component in the implementation of MPC and the impact of modeling uncertainty in MPC is completely contained in this step. Figure 4.10. Comparison of output uncertainty for EnergyPlus simulations with simulated input uncertainty in three input categories, top 10 parameters and top 5 parameters based on the DGSM indices. 66 Other modes of uncertainty that impact MPC operation pertain to real data measurements and feedback, which isn’t the focus of this study. Hence, the brief look into building performance optimization uncertainty that is provided here is used as a platform for a discussion on the impact that uncertainty may have in an MPC framework. The building performance optimization objective, in this case, was to obtain constant setpoints that minimize either HVAC energy demand or PMV. The usage of constant setpoints is a simplification in comparison to the actual setpoint schedule used in the EnergyPlus model. If this method was applied to an MPC framework, the setpoints would be optimized at an hourly interval. This simplification was made here in order to isolate the relationship between uncertainty in our parameter sets and the solution to an arbitrary optimization problem. This relationship is shown in Figure 4.11, where the Pareto frontiers resulting from the multiobjective optimization procedure discussed in Section 4.3 are compared. Five different Pareto frontiers are shown and correspond to the baseline, which has no uncertainty adjustment, and the top 10 and top 5 groups for PMV and FEE from Section 4.4.2, which were optimized with 20% uncertainty in the respective parameter sets. A Pareto frontier represents the set of points that cannot be improved further in one objective without degrading the other objective. And this is the set of points that the optimal control law would select from if it were being used in an MPC algorithm. Four different trials are displayed in Figure 4.11 to show the retest reliability for optimization of each model set. What we can see from this is that the multiobjective optimization isn’t repeatable when uncertainty is applied to the building model inputs, as it is for the baseline case where no uncertainty is applied. This alone could be of concern with regard to the performance of MPC in a real buildings for the same reason that nonconvex problems are a challenge in real-world optimization. Based on the low reliability of the calculation of optimal outputs, the uncertainty contained in the parameters selected here (top 5 and top 10) have made it difficult to obtain a repeatable solution. To confirm this assertion, we would also want to look at the optimal temperature setpoints to identify the difference between solutions for different trials within the same domain. Another interesting observation in Figure 4.11 is that the baseline model acts as an upper bound for the objective function evaluations for all of the optimal solutions. In other words, the optimal performance is always more optimistic when uncertainty is considered. This 67 Figure 4.11. Comparison of Pareto frontiers achieved with the optimization procedure repeated four times for five different uncertainty scenarios: no uncertainty and 20% uncertainty applied to the top 5 and 10 parameters for the PMV and the FEE output. occurs because the spread of the candidate points associated with all of the sets with applied uncertainty will tend to be larger than the spread of the baseline set. And, since the Pareto frontier is a group of extrema from these sets, the Pareto frontiers for the sets that contain uncertainty will always lie underneath the baseline unless the means are significantly different. This could disrupt the connection between the desired and effective weight of objective functions in the selection of an optimal control law in MPC. Since control laws in MPC are calculated based on performance based on a set of objective functions, if the performance of one objective becomes inaccurately optimistic then weight of that objective would effectively be reduced. For example, the PMV Opt 10 Pareto frontier for Trial 1 is almost completely dominated by the thermal comfort objective (i.e., small changes in HVAC extrema give large changes in thermal comfort extrema). This could result in the algorithm putting more effort into minimizing the thermal comfort objective, which doesn’t reflect the true relationship shown in the baseline case. 68 4.5 Summary When performing a UA or SA analysis on a large computational model, it is typically expected that many of the input parameters have a small impact on the model output [101]. This is no exception in this case, where we see that the majority of the parameters considered in the initial uncertainty analysis did not protrude past a noticeable threshold in the sensitivity analysis. The parameters belonging to the materials group had a limited impact on the output uncertainty (Figure 4.5) and only four material parameters made it into the list of high impact parameters shown in Table 4.8 and Table 4.9. Here, we only perturbed the input parameters within an uncertainty interval of 20%, so this by no means indicates that these parameters can be neglected in the model. However, it would not be economical to prioritize reducing the uncertainty in material parameters. Rather, it would be more effective to focus on addressing uncertainty in parameters listed in Table 4.8 or Table 4.9 if attempting to limit uncertainty in FEE or PMV predictions, respectively. It is important to note that the GPR model did not perfectly represent the EnergyPlus model for both output variables. In all of the GPR models there was significant error in the output distribution variance (Table 4.7 and Table 4.11). While this has a smaller impact on derivative-based global sensitivity indices compared to variance based techniques, it certainly still has a nonnegligible impact on the calculation of the DGSM indices [91]. However, the results that were obtained here were consistent and intuitively sensible, so this is not identified as a catastrophic issue, just something that should be addressed in future work. When evaluating the high impact parameters listed in Table 4.8 and Table 4.9, it is notable that the global sensitivity of the rank one parameter for both outputs is on the order of 10 to 100 times greater than the sequential rankings in both cases, which should prompt further investigation into these parameters. In the case of parameter 188 (Activity Schedule) which, numerically, is the metabolic rate in W/m2 , we can attribute its significance directly to its importance in the PMV calculation that is performed by EnergyPlus. For the FEE output, the results are more interesting. In particular, parameter 853 would not have previously been expected to have the highest impact on the output. This parameter is the off cycle parasitic fuel consumption rate (i.e., the electricity consumption of the water heater when it isn’t actively heating water). Most of the remaining parameters listed in Table 4.8 are high 69 energy consumers in the building model, which could be a direct indicator of their influence on output uncertainty in this case. These differences can be explained by the architecture of EnergyPlus and the methodology that it uses to compute different outputs. The methods used here can also indirectly identify known and possibly unknown limitations within the analytic model. Necessarily, this is at the heart of any building modeler’s job and highlights the importance of engineering judgment and continuous scrutiny when working with complex energy system models in any field. By looking at the behavior of optimal solution evaluations when uncertainty in various parameters is present, it is evident that the impact of uncertainty cannot be overlooked. Uncertainty in model inputs can impact the operation of an MPC scheme by convoluting the connection between the optimal inputs and the expected versus actual output performance. The specific challenges that this poses and the associated control behavior will be the subject of future work. CHAPTER 5 CONCLUDING REMARKS Technological advancement, social change, and climate change are all prominent factors that have influenced the evolution of our human experience in the past century. All of these things are interconnected, and in many cases the natural evolution that drives them forward together are irreversible. This is the case today, where simultaneously we have seen tremendous improvements in social standards and quality of life while a large portion of our natural spaces have been altered, contaminated and void of the ecological services we previously relied on. Future efforts will only seek to improve living standards, and as the global population increases, the demand for natural resources that are essential in supporting the modern lifestyle is only expected to grow. So if our aim is to sustain this improvement in living standards into the foreseeable future, then new technology must fill the gap between increasing demand and diminishing resources and a changing climate. There are many sectors of human and physical infrastructure that will need to adapt to changes in climate and resource availability. Agriculture, transportation, energy, water, health, and security are all susceptible to the changes that we expect to see. This thesis focuses on building energy systems, which directly and indirectly impact energy, water, health, and perhaps most importantly, human comfort and productivity. Accounting for up to 40% of global energy consumption, building energy management is now recognized as a field that has significant potential in contributing to climate change and resource management solutions. With this recognition, standards, mandates, and other programs that push forward positive energy management practices are being adopted around the world. For example, many states and cities in the U.S. have adopted benchmarking policies that push building managers to track energy usage and consider ways in which energy performance can be improved [102]. In particular, the University of Utah Energy Management Team has been executing a building metering initiative that will provide 71 access to data that will support the implementation of building control techniques such as MPC in the near future. They plan to use these strategies to meet the University’s carbon reduction goals. MPC in building energy systems has been in a research and development phase for many years, and the key challenges that need to be faced in both building performance simulation (BPS) and MPC itself before it can be implemented on a wide scale are well documented [58], [59], [13]. The complexity associated with modeling a building in a way that supports real-time simulation and control while retaining accuracy proves to be a topical challenge associated with model-based control. Chapter 4 shows us that the modeling approach can be simplified by either using reduced order models or by focusing modeling efforts on areas that have a direct impact on end performance. The utilization of BPS with continuous operation and control will demand that the typical performance gaps in BPS be close to eliminated, while the actual adoption of model-based control requires that the modeling component is technically manageable. There is a dichotomy between the traditional approach to solving these two things separately, making this a problem that demands creative and careful evaluation. As the frequency at which new methods of model-based control are implemented in real buildings increases, the knowledge base related to its use will become an important reference point. So it is essential that researchers and practitioners continue to document case studies as has been done in recent years [103], [104]. The high interest in model-based control techniques is exciting within the context of building simulation and control. But it is also an exciting and necessary component of an interconnected electrical distribution system that can manage fluctuating power supplies from renewable resources. These types of advancements in our energy infrastructure will allow buildings to operate generally as a responsive load, effectively contributing to a more efficient grid infrastructure. In a changing environmental era, where environmental and resource constraints are central to the advancement of new technology and infrastructure, methods that increase the flexibility and responsiveness of the demand base is essential. And, if implemented on a large scale, they can allow us to significantly reduce the total environmental impact of the built environment while supporting a high standard of living for a growing population. REFERENCES [1] WattTime, “WattTime API,” https://api.watttime.org, 2016. [2] J. Macknick, R. Newmark, G. Heath, and K. C. Hallett, “Operational water consumption and withdrawal factors for electricity generating technologies: A review of existing literature,” Environ. Res. Lett., vol. 7, no. 4, pp. 1–10, 2012. “Emissions and generation [3] U.S. Environmental Protection Agency, resource integrated database (eGRID),” https://www.epa.gov/energy/ emissions-generation-resource-integrated-database-egrid, 2014. [4] North American Electric Reliability Corporation, “Key players,” http://www.nerc. com/AboutNERC/keyplayers/Pages/default.aspx, 2017. [5] M. J. Bradley & Associates, “Benchmarking air emissions of the largest electric power producers in the United States,” Tech. Rep., 2017. [6] C. Peng, “Calculation of a building’s life cycle carbon emissions based on Ecotect and building information modeling,” J. Cleaner Prod., vol. 112, pp. 453–465, 2016. [7] M. A. Maupin, J. F. Kenny, S. S. Hutson, J. K. Lovelace, N. L. Barber, and K. S. Linsey, “Estimated use of water in the United States in 2010,” US Geological Survey, Tech. Rep., 2010. [8] D. R. Tobergte and S. Curtis, “DOE Water-Energy Nexus 2014,” J. Chem. Inf. Model., vol. 53, no. 9, pp. 1689–1699, 2013. [9] EIA, “EIA International Energy Outlook 2017,” Energy Information Administration, Tech. Rep., 2017. [10] L. Pérez-Lombard, J. Ortiz, and C. Pout, “A review on buildings energy consumption information,” Energy Build., vol. 40, no. 3, pp. 394–398, 2008. [11] D. B. Crawley, J. W. Hand, M. Kummert, and B. T. Griffith, “Contrasting the capabilities of building energy performance simulation programs,” Build. Environ., vol. 43, no. 4, pp. 661–673, 2008. [12] U.S. Department of Energy, “EnergyPlus,” https://energyplus.net/, 2019. [13] T. Hong, J. Langevin, and K. Sun, “Building simulation: Ten challenges,” Build. Simul., vol. 11, no. 5, pp. 871–898, 2018. [14] S. C. Bengea, A. D. Kelman, F. Borrelli, and R. Taylor, “Model predictive control for mid-size commercial building HVAC: Implementation, results and energy savings,” in Proceedings: The Second International Conference on Building Energy and Environment, 2012, pp. 979–986. 73 [15] M. Morari and J. H. Lee, “Model predictive control: Past, present and future,” Comput. Chem. Eng., vol. 23, pp. 667–682, 1999. [16] D. Mayne, J. Rawlings, C. V. Rao, and P. O. M. Scokaert, “Constrained model predictive control: Stability and optimality,” Automatica, vol. 36, pp. 789–814, 2000. [17] F. Borrelli, Constrained optimal control of linear and hybrid systems, ser. Lecture Notes in Control and Information Sciences. Springer Berlin Heidelberg, 2003. [Online]. Available: https://books.google.com/books?id=mEboCAAAQBAJ [18] R. K. Pachauri, L. Meyer, and The Core Writing Team, Intergovernmental Panel on Climate Change, Tech. Rep. [19] R. Heinberg and D. Lerch, The Post Carbon Reader: Managing the 21st century’s sustainability crises. Watershed Media, 2010. [20] U.S. Environmental Protection Agency, U.S. Environmental Protection Agency, Tech. Rep. [21] J. J. Peck and A. D. Smith, “Modeling power generation water usage,” in Proceedings of the ASME 2015 Power and Energy Conversion Conference, 2015. [22] K. T. Sanders, “Critical review: Uncharted waters? The future of the electricity-water nexus,” Environ. Sci. Technol., vol. 49, no. 1, pp. 51–66, 2015. [23] R. Peer, J. Garrison, C. Timms, and K. Sanders, “Spatially and temporally resolved analysis of environmental trade-offs in electricity generation,” Environ. Sci. Technol., vol. 50, no. 8, pp. 4537–4545, 2016. [24] K. Siler-Evans, I. L. Azevedo, M. G. Morgan, and J. Apt, “Regional variations in the health, environmental, and climate benefits of wind and solar generation.” vol. 110, no. 29. National Academy of Sciences, 2013, pp. 68–73. [25] D. S. Callaway, M. Fowlie, and G. Mccormick, “Location, location, location: The variable value of renewable energy and demand-side efficiency resources,” Environ. Resour. Econ., vol. 5, no. 1, pp. 39–75, 2018. [26] J. S. Graff Zivin, M. J. Kotchen, and E. T. Mansur, “Spatial and temporal heterogeneity of marginal emissions: Implications for electric cars and other electricity-shifting policies,” J. Econ. Behav. Organ., vol. 107, pp. 248–268, 2014. [27] M. M. Rogers, Y. Wang, C. Wang, S. P. McElmurry, and C. J. Miller, “Evaluation of a rapid LMP-based approach for calculating marginal unit emissions,” Appl. Energy, vol. 111, pp. 812–820, 2013. [28] I. Shavel, M. Hagerty, N. Powers, and Y. Yang, “Electricity baseline report for the U.S. power system,” The Brattle Group, Tech. Rep., 2015. [29] K. Siler-Evans, I. L. Azevedo, and M. G. Morgan, “Marginal emissions factors for the U.S. electricity system,” Environ. Sci. Technol., vol. 46, no. 9, pp. 4742–4748, 2012. [30] J. Lariviere, G. Mccormick, S. M. Browder, S. M. Harriman, M. Murray, and S. Keeton, “How better accounting can more cheaply reduce carbon emissions,” The Howard H. Baker Jr. Center for Public Policy, Tech. Rep., 2016. 74 [31] M. Sotos, “GHG protocol scope 2 guidance: An amendment to the GHG protocol corporate standard,” World Resources Institute, Tech. Rep., 2015. [32] U.S. Environmental Protection Agency, “Air Markets Program Data–Clean Air Markets–US EPA,” https://ampd.epa.gov/ampd, 2018. [33] G. Pétron, P. Tans, G. Frost, D. Chao, and M. Trainer, “High-resolution emissions of CO2 from power generation in the USA,” J. Geophys. Res. Biogeosci., vol. 113, no. 4, pp. 1–9, 2008. [34] J. A. Gouw, D. D. Parrish, G. J. Frost, and M. Trainer, “Reduced emissions of CO2, NOx, and SO2 from U.S. power plants owing to switch from coal to natural gas with combined cycle technology,” Earth’s Future, vol. 2, no. 2, pp. 75–82, 2014. [35] U.S. Energy Information Administration, “EIA Open Data–Bulk–U.S. Energy Information Administration,” https://www.eia.gov/opendata/bulkfiles.php, 2018. [36] G. Schivley, I. Azevedo, and C. Samaras, “Assessing the evolution of power sector carbon intensity in the United States,” Environ. Res. Lett., vol. 13, no. 6, 2018. [37] J. F. Kenny, N. L. Barber, S. S. Hutson, K. S. Linsey, J. K. Lovelace, and M. A. Maupin, “Estimated use of water in the United States in 2005,” US Geological Survey, Tech. Rep., 2009. [38] B. K. Sovacool and K. E. Sovacool, “Identifying future electricity-water tradeoffs in the United States,” Energy Policy, vol. 37, no. 7, pp. 2763–2773, 2009. [39] F. Ackerman and J. Fisher, “Is there a water-energy nexus in electricity generation? Long-term scenarios for the western United States,” Energy Policy, vol. 59, pp. 235–241, 2013. [40] K. M. Strzepek, J. Baker, W. Farmer, and C. A. Schlosser, “Modeling water withdrawal and consumption for electricity generation in the United States,” MIT Joint Program on the Science and Policy of Global Change, Tech. Rep. 222. [41] J. Macknick, S. Sattler, K. Averyt, S. Clemmer, and J. Rogers, “The water implications of generating electricity: Water use across the United States based on different electricity pathways through 2050,” Environ. Res. Lett., vol. 7, no. 4, pp. 1–10, 2012. [42] J. Macknick, J. Meldrum, S. Nettles-Anderson, G. Heath, and A. Miara, “Life cycle water use for photovoltaic electricity generation: A review and harmonization of literature estimates,” Environ. Res. Lett., vol. 8, pp. 1–18, 2013. [43] N. Bracken, “Water and energy in the west: The legal and institutional issues that affect water availability for energy-related activities,” Western States Water Council, Tech. Rep., 2015. [44] C. A. Scott, S. A. Pierce, M. J. Pasqualetti, A. L. Jones, B. E. Montz, and J. H. Hoover, “Policy and institutional dimensions of the water-energy nexus,” Energy Policy, vol. 39, no. 10, pp. 6622–6630, 2011. [45] U.S. Environmental Protection Agency, “Power Profiler,” https://www.epa.gov/ energy/power-profiler, 2017. 75 [46] U.S. Energy Information Administration, “Electricity Data Browser,” https://www. eia.gov/electricity/data/browser, 2017. [47] Site-Specific Energy Systems Laboratory Research Group, “Electricity and Environment Repository,” 2018. [Online]. Available: https://github.com/SSESLab/ electricityandenvironment [48] S. Clemmer, J. Rogers, S. Sattler, J. Macknick, and T. Mai, “Modeling low-carbon US electricity futures to explore impacts on national and regional water use,” Environ. Res. Lett., vol. 8, no. 1, pp. 1–11, 2013. [49] K. Janda, “Worldwide status of energy standards for buildings,” Energy, vol. 19, no. 1, pp. 27–44, 2009. [50] T. I. Salsbury, “A survey of control technologies in the building automation industry,” in Proceedings of the 16th Triennial World Congress, vol. 16. IFAC, 2005, pp. 90–100. [51] V. S. Harish and A. Kumar, “A review on modeling and simulation of building energy systems,” Renewable Sustainable Energy Rev., vol. 56, pp. 1272–1292, 2016. [52] E. Shen and T. Hong, “Simulation-based assessment of the energy savings benefits of integrated control in office buildings,” Build. Simul., vol. 2, no. 4, pp. 239–251, 2009. [53] F. Oldewurtel, A. Parisio, C. N. Jones, D. Gyalistras, M. Gwerder, V. Stauch, B. Lehmann, and M. Morari, “Use of model predictive control and weather forecasts for energy efficient building climate control,” Energy Build., vol. 45, pp. 15–27, 2012. [54] G. P. Henze and M. Krarti, “Predictive optimal control of active and passive building thermal storage inventory,” HVAC&R Res., vol. 11, no. 2, pp. 189–213, 2005. [55] S. H. Cho and M. Zaheer-Uddin, “Predictive control of intermittently operated radiant floor heating systems,” Energy Convers. Manage., vol. 44, no. 8, pp. 1333–1342, 2003. [56] D. Bourgeois, C. Reinhart, and I. Macdonald, “Adding advanced behavioural models in whole building energy simulation: A study on the total energy impact of manual and automated lighting control,” Energy Build., vol. 38, no. 7, pp. 814–823, 2006. [57] F. Oldewurtel, D. Gyalistras, M. Gwerder, C. Jones, A. Parisio, V. Stauch, B. Lehmann, and M. Morari, “Increasing energy efficiency in building climate control using weather forecasts and model predictive control,” Clima - RHEVA World Congress, 2010. [58] M. Killian and M. Kozek, “Ten questions concerning model predictive control for energy efficient buildings,” Build. Environ., vol. 105, pp. 403–412, 2016. [59] J. Cigler, D. Gyalistras, J. Siroky, V.-N. Tiet, and L. Ferkl, “Beyond theory: The challenge of implementing model predictive control in buildings,” Clima 2013 - 11th REHVA World Congress & 8th International Conference on IAQVEC - Energy Efficient, Smart and Healthy Buildings, pp. 1008–1018, 2013. [60] U.S. Department of Energy, “Commercial Prototype Building Models,” https://www. energycodes.gov/development/commercial/prototype_models, 2018. 76 [61] M. Liu and Y. Shi, “Model predictive control for thermostatically controlled appliances providing balancing service,” IEEE Trans. Control Syst. Technol., vol. 24, no. 6, pp. 2082–2093, 2016. [62] F. Belić, Ž. Hocenski, and D. Slišković, “Thermal modeling of buildings with RC method and parameter estimation,” in Proceedings of 2016 International Conference on Smart Systems and Technologies, SST 2016. IEEE, 2016, pp. 19–25. [63] F. Wang, “Model predictive control of thermal dynamics in building,” Master’s Thesis, Lehigh University, Bethlehem, PA 18015, 2012. [64] I. Hazyuk, C. Ghiaus, and D. Penhouet, “Optimal temperature control of intermittently heated buildings using model predictive control: Part I - building modeling,” Build. Environ., vol. 51, pp. 379–387, 2012. [65] M. D. Knudsen, R. E. Hedegaard, T. H. Pedersen, and S. Petersen, “System identification of thermal building models for demand response - A practical approach,” Energy Procedia, vol. 122, pp. 937–942, 2017. [66] S. Wilcox and W. Marion, “Users manual for TMY3 data sets,” National Renewable Energy Laboratory, Tech. Rep., 2008. [67] M. Grant and S. Boyd, “CVX: Matlab software for disciplined convex programming, version 2.1,” http://cvxr.com/cvx, 2014. [68] E. Vrettos, J. L. Mathieu, and G. Andersson, “Control of thermostatic loads using moving horizon estimation of individual load states,” in Proceedings - 2014 Power Systems Computation Conference, PSCC 2014. IEEE, 2014. [69] J. Clarke and T. Maver, “Advanced design tools for energy conscious building design: Development and dissemination,” Build. Environ., vol. 26, no. 1, pp. 25 – 34, 1991. [70] I. Macdonald, “Quantifying the effects of uncertainty in building simulation,” PhD Dissertation, University of Strathclyde, 2002. [71] D. Coakley, “Calibration of detailed building energy simulation models to measured data using uncertainty analysis,” PhD Dissertation, NUI Galway, 2014. [72] D. Bronson, S. Hinchey, J. Haberl, and D. O’Neal, “Procedure for calibrating the DOE-2 simulation program to non-weather-dependent measured loads.” in ASHRAE Winter Meeting, 1992, pp. 636–652. [73] F. Campolongo, S. Tarantola, and A. Saltelli, “Sensitivity analysis as an ingredient of modeling,” Stat. Sci., vol. 15, no. 4, pp. 377–395, 2002. [74] B. Eisenhower, Z. O’Neill, S. Narayanan, V. A. Fonoberov, and I. Mezić, “A methodology for meta-model based optimization in building energy models,” Energy Build., vol. 47, pp. 292–301, 2012. [75] B. Eisenhower, O. Zheng, and I. Mezić, “Leveraging the analysis of parametric uncertainty for building energy model calibration,” Build. Simul., vol. 6, no. 4, pp. 365–377, 2013. 77 [76] U.S. Department of Energy, “Input Output Reference,” Tech. Rep., 2019. [77] Y. H. Yau and B. T. Chew, “A review on predicted mean vote and adaptive thermal comfort models,” Build. Serv. Eng. Res. Technol., vol. 35, no. 1, pp. 23–35, 2014. [78] U.S. Department of Energy, “Engineering Reference – Occupant Thermal Comfort,” https://bigladdersoftware.com/epx/docs/8-8/engineering-reference/ occupant-thermal-comfort.html, 2019. [79] W. Tian, “A review of sensitivity analysis methods in building energy analysis,” Renewable Sustainable Energy Rev., vol. 20, pp. 411–419, 2013. [80] T. Østergård, R. L. Jensen, and S. E. Maagaard, “A comparison of six metamodeling techniques applied to building performance simulations,” Appl. Energy, vol. 211, pp. 89–103, 2018. [81] A. Rahman and A. D. Smith, “Predicting fuel consumption for commercial buildings with machine learning algorithms,” Energy and Build., vol. 152, pp. 341–358, 2017. [82] Y. Heo and V. M. Zavala, “Gaussian process modeling for measurement and verification of building energy savings,” Energy and Build., vol. 53, pp. 7–18, 2012. [83] B. Eisenhower, Z. O’Neill, V. A. Fonoberov, and I. Mezić, “Uncertainty and sensitivity decomposition of building energy models,” J. Build. Perform. Simul., vol. 00, pp. 1–18, 2011. [84] M. Manfren, N. Aste, and R. Moshksar, “Calibration and uncertainty analysis for computer models - A meta-model based approach for integrated building energy simulation,” Appl. Energy, vol. 103, pp. 627–641, 2013. [85] S. Paudel, P. H. Nguyen, W. L. Kling, M. Elmitri, B. Lacarrière, and O. L. Corre, “Support vector machine in prediction of building energy demand using pseudo dynamic approach,” in Proceedings of ECOS 2015-The 28th International Conference on Efficiency, Cost, Optimization, Simulation and Environmental Impact of Energy Systems. ECOS, 2015. [86] Y. Zhang, Z. O’Neill, B. Dong, and G. Augenbroe, “Comparisons of inverse modeling approaches for predicting building energy performance,” Build. Environ., vol. 86, pp. 177–190, 2015. [87] F. Pedregosa, G. Varoquaux et al., “Scikit-learn: Machine learning in Python,” https: //scikit-learn.org/stable/index.html, 2007. [88] F. Pedregosa and G. Varoquaux, “Scikit-learn: Machine learning in Python,” J. Mach. Learn. Res., vol. 12, pp. 2825–2830, 2011. [89] R. Venâncio and A. Pedrini, “The influence of design decisions on energy consumption and thermal performance: The case of UFRN campus, Brazil,” Build. Simul., pp. 137–143, 2009. [90] M. Pollock, Y. Roderick, D. Mcewan, and C. Wheatley, “Building simulation as an assisting tool in designing an energy efficient building: A case study,” Build. Simul., pp. 1191–1198, 2009. 78 [91] I. M. Sobol’ and S. Kucherenko, “Derivative based global sensitivity measures and their link with global sensitivity indices,” Math. Comput. Simul., vol. 79, no. 10, pp. 3009–3017, 2009. [92] J. Herman, W. Usher et al., “SALib: Sensitivity analysis library in Python,” https: //salib.readthedocs.io/en/latest/index.html, 2017. [93] C. E. Rasmussen and C. K. I. Williams, Gaussian processes for machine learning. MIT Press, 2006. [Online]. Available: www.GaussianProcess.org/gpml The [94] M. Lamboni, B. Iooss, A. L. Popelin, and F. Gamboa, “Derivative-based global sensitivity measures: General links with Sobol’ indices and numerical tests,” Math. Comput. Simul., vol. 87, pp. 45–54, 2013. [95] D. Hadka, “Platypus Documentation Contents,” https://buildmedia.readthedocs. org/media/pdf/platypus/latest/platypus.pdf, 2018. [96] K. Deb, A. Pratap, S. Agarwal, and T. Meyarivan, “A fast and elitist multi-objective genetic algoritm:NSGA -II,” IEEE Trans. Evol. Comput., vol. 6, no. 2, pp. 182–197, 2002. [97] S. Philip, “eppy Documentation,” https://buildmedia.readthedocs.org/media/pdf/ eppy/latest/eppy.pdf, 2018. [98] K. Plewe and A. D. Smith, “PECT: A tool for computing the temporal and spatial variation of externalities related to power generation in the United States,” SoftwareX, vol. 9, pp. 61–67, 2019. [99] E. Local, “Electricity Local: Local electricity rates and statistics,” https://www. electricitylocal.com/, 2019. [100] A. Bhattacharyya, “On a measure of divergence between two statistical populations defined by their probability distributions,” Bull. Calc. Math. Soc., vol. 35, pp. 99–110, 1943. [101] A. Saltelli, K. Chan, M. Scott et al., Sensitivity analysis. Probability and statistics series. New York: John Wiley & Sons, 2000. [102] U.S. Environmental Protection Agency, “Benchmarking programs and policies leveraging Energy Star benchmarking programs and policies,” Tech. Rep., 2019. [103] A. Afram, F. Janabi-Sharifi, A. S. Fung, and K. Raahemifar, “Artificial neural network (ANN) based model predictive control (MPC) and optimization of HVAC systems: A state of the art review and case study of a residential HVAC system,” Energy Build., vol. 141, pp. 96–113, 2017. [104] D. Sturzenegger, D. Gyalistras, M. Morari, and R. S. Smith, “Model predictive climate control of a Swiss office building,” IEEE Trans. Control Syst. Technol., vol. 24, no. 1, pp. 1–12, 2016.