Patent Publication Number: US-10325331-B2

Title: Systems and methods for measuring and verifying energy usage in a building

Description:
CROSS-REFERENCE TO RELATED PATENT APPLICATION 
     This application is a continuation of U.S. patent application Ser. No. 13/485,682 filed May 31, 2012, the entire disclosure of which is incorporated by reference herein. 
    
    
     BACKGROUND 
     The present disclosure generally relates to energy conservation in a building. The present disclosure relates more specifically to detecting changes in the energy use of a building. 
     In many areas of the country, electrical generation and transmission assets have or are reaching full capacity. One of the most cost effective ways to ensure reliable power delivery is to reduce demand (MW) by reducing energy consumption (MWh). Because commercial buildings consume a good portion of the generated electricity in the United States, a major strategy for solving energy grid problems is to implement energy conservation measures (ECMs) within buildings. Further, companies that purchase energy are working to reduce their energy costs by implementing ECMs within buildings. 
     Entities that invest in ECMs typically want to verify that the expected energy savings associated with ECMs are actually realized (e.g., for verifying the accuracy of return-on-investment calculations). Federal, state, or utility based incentives may also be offered to encourage implementation of ECMs. These programs will have verification requirements. Further, some contracts between ECM providers (e.g., a company selling energy-efficient equipment) and ECM purchasers (e.g., a business seeking lower ongoing energy costs) establish relationships whereby the ECM provider is financially responsible for energy or cost savings shortfalls after purchase of the ECM. Accordingly, Applicants have identified a need for systems and methods for measuring and verifying energy savings and peak demand reductions in buildings. Applicants have further identified a need for systems and methods that automatically measure and verify energy savings and peak demand reductions in buildings. 
     Once a baseline model has been developed for measuring and verifying energy savings in a building, the baseline model may be used to estimate the amount of energy that would have been used by the building if ECMs were not performed. Calculation of the energy avoidance may be rather simple if the baseline model is accurate and routine adjustments related to building usage are included in the baseline model. However, non-routine adjustments related to building usage may not be included in the baseline model, and therefore not included in the energy avoidance calculation. Examples of static factors for which non-routine adjustments may occur include the size of the facility, the hours of operation of the building, the number of employees using the building, the number of computer servers in the building, etc. 
     Such non-routine adjustments may cause problems during the reporting period when trying to verify energy savings in a building. An unnoticed change (e.g., a non-routing adjustment) in static factors may cause energy savings in the building to be underestimated. This may cause an engineer to miss requirements in a contract or cause a building owner not to invest money in ECMs even though the actual return on investment was high. 
     Unaccounted for changes in static factors that occur during the baseline period may also negatively impact calculations (e.g., creation of an accurate baseline model). If changes in static factors occur during the baseline period, the statistical regression module used to build the baseline model may essentially be based on a midpoint between energy usage before the change in static factors and energy usage after the change in static factors. Such a situation may lead to an inaccurate estimation of energy usage. For example, if a change in static factors causes a decrease in energy usage, the energy savings of the building due to the ECMs will be overestimated by the baseline model, and the building owner may be disappointed with a resulting return on investment. As another example, if a change in static factors causes an increase in energy usage, the energy savings of the building will be underestimated by the baseline model, and a performance contractor may miss the contractual guarantee of energy savings and be required to reimburse the building owner. Monitoring static factors traditionally requires manual monitoring. This may lead to changes in static factors going unnoticed or incorrectly calculated. 
     SUMMARY 
     One embodiment of the invention relates to a computerized method of determining the end of a measurement and verification period of a building management system. The method includes calculating, for a plurality of periods, a target amount of energy savings resulting from energy conservation measures in a building. The method further includes calculating, for a plurality of periods, a least amount of energy savings resulting from energy conservation measures in a building. The method further includes comparing the target amount to the least amount. The method further includes ending the measurement and verification period, before the end of a maximum measurement and verification period, when the least amount is greater than the target amount. The method further includes outputting an indication that the measurement and verification period has ended to at least one of a memory device, a user device, or another device on the building management system. 
     Another embodiment of the invention relates to a computerized method of determining the end of a measurement and verification period of a building management system. The method includes calculating, for a plurality of periods, a target amount of energy savings resulting from energy conservation measures in a building. The method further includes calculating, for a plurality of periods, a greatest amount of energy savings resulting from energy conservation measures in a building. The method further includes comparing the target amount to the greatest amount. The method further includes ending the measurement and verification period, before the end of a maximum measurement and verification period, when the greatest amount is less than the target amount. The method further includes outputting an indication that the measurement and verification period has ended to at least one of a memory device, a user device, or another device on the building management system. 
     Yet another embodiment of the invention relates to a controller for determining the end of a measurement and verification period of a building management system. The controller includes a processing circuit configured to calculate, for a plurality of periods, a target amount of energy savings resulting from energy conservation measures in a building. The processing circuit is further configured to calculate, for a plurality of periods, a least amount of energy savings resulting from energy conservation measures in a building. The processing circuit is further configured to compare the target amount to the least amount. The processing circuit is further configured to end the measurement and verification period, before the end of a maximum measurement and verification period, when the least amount is greater than the target amount. The processing circuit is further configured to output an indication that the measurement and verification period has ended to at least one of a memory device, a user device, or another device on the building management system. 
     Yet another embodiment of the invention relates to computer-readable media with computer-executable instructions embodied thereon that when executed by a computer system perform a method for use with a building management system in a building. The instructions include instructions for calculating, for a plurality of periods, a target amount of energy savings resulting from energy conservation measures in a building. The instructions further include instructions for calculating, for a plurality of periods, a least amount of energy savings resulting from energy conservation measures in a building. The instructions further include instructions for comparing the target amount to the least amount. The instructions further include instructions for ending the measurement and verification period, before the end of a maximum measurement and verification period, when the least amount is greater than the target amount. The instructions further include instructions for outputting an indication that the measurement and verification period has ended to at least one of a memory device, a user device, or another device on the building management system. 
     Alternative exemplary embodiments relate to other features and combinations of features as may be generally recited in the claims. 
    
    
     
       BRIEF DESCRIPTION OF THE FIGURES 
       The disclosure will become more fully understood from the following detailed description, taken in conjunction with the accompanying figures, wherein like reference numerals refer to like elements, in which: 
         FIG. 1A  is a flow chart of a process for measuring and verifying energy savings and peak demand reductions in a building, according to an exemplary embodiment; 
         FIG. 1B  is a simplified block diagram of a system for completing or facilitating the process of  FIG. 1A , according to an exemplary embodiment; 
         FIG. 10  is a block diagram of a system for measuring and verifying energy savings in a building, according to an exemplary embodiment; 
         FIG. 2  is a detailed block diagram of the baseline calculation module of  FIG. 10 , according to an exemplary embodiment; 
         FIG. 3A  is a flow chart of a process for selecting observed variable data to use for generation of the baseline model, according to an exemplary embodiment; 
         FIG. 3B  is a flow chart of a process for selecting calculated variable data to use for generation of the baseline model, according to an exemplary embodiment; 
         FIGS. 4A-4E  are more detailed flow charts of the process of  FIG. 3B , according to an exemplary embodiment; 
         FIG. 5  is a flow chart of the objective function used in the golden section search of the process of  FIGS. 4A-E  shown in greater detail, according to an exemplary embodiment; 
         FIG. 6  is a flow chart of a process of calculating enthalpy, according to an exemplary embodiment; 
         FIG. 7  is a flow chart of a process for identifying changes in the static factors of a baseline model such as that discussed with reference to previous Figures, according to an exemplary embodiment; 
         FIG. 8A  is a flow chart of a process for detecting changes in static factors of a baseline model during a baseline period, according to an exemplary embodiment; 
         FIG. 8B  is a flow chart of a process for calculating an observed statistic for detecting changes in static factors during a baseline period, according to an exemplary embodiment; 
         FIG. 9A  is a flow chart of a process for detecting changes in static factors of a baseline model during a reporting period, according to an exemplary embodiment; 
         FIG. 9B  is a flow chart of a process for calculating an observed statistic for detecting changes in static factors during a reporting period, according to an exemplary embodiment; 
         FIG. 10  is a detailed block diagram of a baseline calculation module for detecting changes during a baseline period, according to an exemplary embodiment; 
         FIG. 11  is a detailed block diagram of a facility monitoring module for detecting changes during a reporting period, according to an exemplary embodiment; 
         FIG. 12  is a flow chart of evaluating an inverse of a cumulative probability distribution function for determining a critical value for the process of  FIG. 9A , according to an exemplary embodiment; 
         FIG. 13  is a flow chart of calculating a P-value for the process of  FIG. 8A , the P-value for use in determining changes in static factors, according to an exemplary embodiment; 
         FIG. 14  is a block diagram of a system for using uncertainty to shorten the measurement and verification process, according to an exemplary embodiment; 
         FIG. 15  is a detailed block diagram of the baseline calculation module of  FIG. 14 , according to an exemplary embodiment; 
         FIG. 16  is a flow chart for using uncertainty to shorten the measurement and verification process, according to an exemplary embodiment; 
         FIG. 17  is a more detailed flow chart of the process of  FIG. 16 , according to an exemplary embodiment; 
         FIG. 18  is a chart of energy savings and related quantities, according to an exemplary embodiment; 
         FIG. 19A  is a plot of energy savings, lower confidence bounds, upper confidence bounds, and uncertainty as a percentage of energy savings, according to an exemplary embodiment; 
         FIG. 19B  is a plot of target energy savings and lower confidence bounds, according to an exemplary embodiment; and 
         FIG. 19C  is a plot of target energy savings and actual energy savings, according to an exemplary embodiment. 
     
    
    
     DESCRIPTION 
     Before turning to the figures, which illustrate the exemplary embodiments in detail, it should be understood that the disclosure is not limited to the details or methodology set forth in the description or illustrated in the figures. It should also be understood that the terminology is for the purpose of description only and should not be regarded as limiting. 
     Referring generally to the Figures, a computer system for use with a building management system in a building is provided for automatically detecting changes in factors correlated to energy use that are expected to remain static (e.g., one or more static factors). The computer system includes a processing circuit configured to automatically identify a change in such a static factor of a building based on data received from the building management system. The processing circuit may be configured to communicate the identified change in the static factor to at least one of (a) a module for alerting a user to changes and (b) a module for initiating adjustment to an energy model for a building. 
     Embodiments of the present disclosure are configured to automatically (e.g., via a computerized process) calculate a baseline model (i.e., a predictive model) for use in measuring and verifying energy savings and peak demand reductions attributed to the implementation of energy conservation measures in building. The calculation of the baseline model may occur by applying a partial least squares regression (PLSR) method to data from a building management system (BMS). The baseline model is used to predict energy consumption and peak demand reductions in a building if an ECM were not installed or used in the building. Actual energy consumption using the ECM is subtracted from the predicted energy consumption to obtain an energy savings estimate or peak demand estimate. 
     The computerized process can utilize many collinear or highly correlated data points from the BMS to calculate the baseline model using the PLSR algorithm. Data clean-up, data synchronization, and regression analysis activities of the computerized process can be used to prepare the data and to tune the baseline model for improved performance relative to the pertinent data. Further, baseline contractual agreements and violations based on the generated baseline model may be determined with a predetermined statistical confidence. 
     To provide improved performance over conventional approaches to baseline calculations, an exemplary embodiment includes the following features: one or more computerized modules for automatically identifying which predictor variables (e.g., controllable, uncontrollable, etc.) are the most important to an energy consumption prediction and a computerized module configured to automatically determine the baseline model using the identified predictor variables determined to be the most important to the energy consumption prediction. 
     While the embodiments shown in the figures mostly relate to measuring and verifying energy consumption and energy savings in a building with the use of expected values as inputs, it should be understood that the systems and methods below may be used to measure and verify peak demand consumption and savings with the use of maximum values as inputs. 
     Measuring and Verifying Energy Savings in Buildings 
     Referring now to  FIGS. 1A and 1B , a process  100  for measuring and verifying energy savings and peak demand in a building is shown, according to an exemplary embodiment. Process  100  is shown to include retrieving historical building and building environment data  120  from a pre-retrofit period (step  102 ). Input variables retrieved in step  102  and used in subsequent steps may include both controllable variables (i.e., variables that may be controlled by a user such as occupancy of an area and space usage) and uncontrollable variables (e.g., outdoor temperature, solar intensity and duration, humidity, other weather occurrences, etc.). 
     Process  100  further includes using the data obtained in step  102  to calculate and select a set of variables significant to energy usage in the building (step  104 ). Step  104  may include calculating variables that may be used to determine energy usage in the building. For example, calculated variables such as cooling degree days, heating degree days, cooling energy days, or heating energy days that are representative of energy usage in the building relating to an outside air temperature and humidity may be calculated. Energy days (cooling energy days and heating energy days) are herein defined as a predictor variable that combines both outside air temperature and outside air humidity. Energy days differ from degree days at least in that the underlying integration is based on the calculated outside air enthalpy. Step  104  may include the selection of a set of calculated variables, variables based on a data set from the data received in step  102 , or a combination of both. For example, the set of variables may include variables associated with a data set of the building (e.g., occupancy and space usage of an area, outdoor air temperature, humidity, solar intensity) and calculated variables (e.g., occupancy hours, degree days, energy days, etc.). Variables and data that are not significant (e.g., that do not have an impact on energy usage in the building) may be discarded or ignored by process  100 . 
     The set of variables is then used to create a baseline model  126  that allows energy usage or power consumption to be predicted (step  106 ). With reference to the block diagram of  FIG. 1B , baseline model  126  may be calculated using a baseline model generator  122  (e.g., a computerized implementation of a PLSR algorithm). 
     Process  100  further includes storing agreed-upon ranges of controllable input variables and other agreement terms in memory (step  108 ). These stored and agreed-upon ranges or terms are used as baseline model assumptions in some embodiments. In other embodiments the baseline model or a resultant contract outcome may be shifted or changed when agreed-upon terms are not met. 
     Process  100  further includes conducting an energy efficient retrofit of building equipment (step  110 ). The energy efficient retrofit may include any one or more process or equipment changes or upgrades expected to result in reduced energy consumption by a building. For example, an energy efficient air handling unit having a self-optimizing controller may be installed in a building in place of a legacy air handling unit with a conventional controller. 
     Once the energy efficient retrofit is installed, process  100  begins obtaining measured energy consumption  130  for the building (step  112 ). The post-retrofit energy consumption  130  may be measured by a utility provider (e.g., power company), a system or device configured to calculate energy expended by the building HVAC system, or otherwise. 
     Process  100  further includes applying actual input variables  124  of the post-retrofit period to the previously created baseline model  126  to predict energy usage of the old system during the post-retrofit period (step  114 ). This step results in obtaining a baseline energy consumption  128  (e.g., in kWh) against which actual energy consumption  130  from the retrofit can be compared. 
     In an exemplary embodiment of process  100 , estimated baseline energy consumption  128  is compared to measured energy consumption  130  by subtracting measured energy consumption  130  during the post-retrofit period from estimated baseline energy consumption  128  (step  116 ). This subtraction will yield the energy savings  132  resulting from the retrofit. The energy savings  132  resulting from the retrofit is multiplied or otherwise applied to utility rate information for the retrofit period to monetize the savings (step  118 ). Steps  114  and  116  may further include determining a peak demand reduction in the building and monetizing cost related to the reduction. 
     Referring now to  FIG. 10 , a more detailed block diagram of a BMS computer system  200  for measuring and verifying energy savings in a building is shown, according to an exemplary embodiment. System  200  includes multiple inputs  202  from disparate BMS sources. Inputs  202  are received and parsed or otherwise negotiated by an information aggregator  204  of the processing circuit  254 . 
     BMS computer system  200  includes a processing circuit  250  including a processor  252  and memory  254 . Processor  252  can be implemented as a general purpose processor, an application specific integrated circuit (ASIC), one or more field programmable gate arrays (FPGAs), a group of processing components, or other suitable electronic processing components. Memory  254  is one or more devices (e.g., RAM, ROM, Flash memory, hard disk storage, etc.) for storing data and/or computer code for completing and/or facilitating the various processes, layers, and modules described in the present disclosure. Memory  254  may be or include volatile memory or non-volatile memory. Memory  254  may include database components, object code components, script components, or any other type of information structure for supporting the various activities and information structures described in the present disclosure. According to an exemplary embodiment, memory  254  is communicably connected to processor  252  via processing circuit  250  and includes computer code for executing (e.g., by processing circuit  250  and/or processor  252 ) one or more processes described herein. 
     Memory  254  includes information aggregator  204 . Information aggregator  204  may serve as middleware configured to normalize communications or data received from the multiple inputs. Information aggregator  204  may be a middleware appliance or module sold by Johnson Controls, Inc. Information aggregator  204  is configured to add data to a BMS database  206 . A data retriever  208  is configured to retrieve (e.g., query) data from BMS database  206  and to pass the retrieved data to baseline calculation module  210 . 
     Baseline calculation module  210  is configured to create a baseline model using historical data from the multiple inputs and aggregated in BMS database  206  or other data sources. Some of the information may be received from sources other than building data sources (e.g., weather databases, utility company databases, etc.). The accuracy of the baseline model will be dependent upon errors in the data received. 
     Baseline calculation module  210  is shown to include data clean-up module  212 . Data clean-up module  212  receives data from data retriever  208  and prefilters the data (e.g., data scrubbing) to discard or format bad data. Data clean-up module  212  conducts one or more checks to determine whether the data is reliable, whether the data is in the correct format, whether the data is or includes a statistical outlier, whether the data is distorted or “not a number” (NaN), whether the sensor or communication channel for a set of data has become stuck at some value, and if the data should be discarded. Data clean-up module  212  may be configured to detect errors via, for example, threshold checks or cluster analysis. 
     Baseline calculation module  210  is further shown to include data synchronization module  214 . Data synchronization module  214  receives the data after the data is “cleaned up” by data clean-up module  212  and is configured to determine a set of variables for use in generating the baseline model. The variables may be calculated by module  214 , may be based on received data from data retriever  208  and data clean-up module  212 , or a combination of both. For example, data synchronization module  214  may determine variables (e.g., cooling and heating degree days and cooling and heating energy days) that serve as a proxy for energy usage needed to heat or cool an area of the building. Data synchronization module  214  may then further determine which type of calculated variable to use (e.g., whether to use degree days or energy days in the regression analysis that generates the baseline model). Further, data synchronization module  214  may identify and use measured data received from data retriever  208  and formatted by data clean-up module for use in the set of variables. For example, module  214  may select temperature data received from data retriever  208  as a predictor variable for energy usage in the building. 
     Baseline calculation module  210  further includes regression analysis module  216 . Regression analysis module  216  is configured to generate the baseline model based on the set of variables from data synchronization module  214 . According to one exemplary embodiment, a partial least squares regression (PLSR) method may be used to generate the baseline model. According to other embodiments, other regression methods (e.g., a principal component regression (PCR), ridge regression (RR), ordinary least squares regression (OLSR)) are also or alternatively used in the baseline model calculation. The PLSR method is based on a linear transformation from the set of variables from module  214  to a linear model that is optimized in terms of predictivity. 
     Baseline calculation module  210  further includes cross-validation module  218 . Cross-validation module  218  is configured to validate the baseline model generated by regression analysis module  216 . Validation of the baseline model may include ensuring there is no overfitting of the baseline model (e.g., having too many variables or inputs influencing the model), determining a correct order or number of components in the model, or conducting other tests or checks of the baseline module output by regression analysis module  216 . Baseline calculation module  210  and sub-modules  212 - 218  are shown in greater detail in  FIG. 2  and subsequent figures. 
     Post retrofit variable data  222  is applied to baseline model  220  generated by baseline calculation module  210  (e.g., data relating to estimated energy use of the building) to obtain a resultant baseline energy consumption. Measured energy consumption  224  from the building is subtracted from the resultant baseline energy consumption at element  225  to obtain energy savings data  226 . Energy savings data  226  may be used to determine payments (e.g., from the retrofit purchaser to the retrofit seller), to demonstrate the new equipment&#39;s compliance with a guaranteed level of performance, or as part of a demand-response commitment or bid validation. Energy savings data  226  may relate to overall energy consumption and/or peak demand reduction. 
     Energy savings data  226  or other information from the prior calculations of the system is used to monitor a building after retrofit to ensure that the facility occupants have not violated the terms of the baseline agreement (e.g., by substantially adding occupants, by changing the building space use, by bringing in more energy using devices, by substantially changing a setpoint or other control setting, etc.). Conventionally this involves making periodic visits to the facilities, reviewing job data, and/or making specialized measurements. Because visits and specialized reviews are time consuming, they are often not done, which puts any savings calculations for a period of time in question. 
     System  200  includes an Exponentially Weighted Moving Average (EWMA) control module  228  configured to automate a baseline term validation process. EWMA control module  228  monitors the difference between the predicted and measured consumption. Specifically, EWMA control module  228  checks for differences between the predicted and measured consumption that are outside of predetermined statistical probability thresholds and provides the results to a facility monitoring module  230 . Any statistically unlikely occurrences can cause a check of related data against baseline agreement information  232 , used to update baseline agreement information, or are provided to a user output/system feedback module  234 . User output/system feedback module  234  may communicate alarm events to a user in the form of a displayed (e.g., on an electronic display) EWMA chart configured to highlight unexpected shifts in the calculated energy savings. Input calibration module  236  may receive feedback from module  234  and additionally provide data to data retriever  208  regarding instructions to add or remove variables from consideration for the baseline model in the future. In other embodiments, different or additional control modules may implement or include statistical process control approaches other than or in addition to EWMA to provide baseline validation features. 
     BMS computer system  200  further includes a user input/output (I/O)  240 . User I/O  240  is configured to receive a user input relating to the data set used by baseline calculation module  210  and other modules of system  200 . For example, user I/O  240  may allow a user to input data into system  200 , edit existing data points, etc. System  200  further includes communications interface  242 . Communications interface  242  can be or include wired or wireless interfaces (e.g., jacks, antennas, transmitters, receivers, transceivers, wire terminals, etc.) for conducting data communications with the BMS, subsystems of the BMS, or other external sources via a direct connection or a network connection (e.g., an Internet connection, a LAN, WAN, or WLAN connection, etc.). 
     Referring now to  FIG. 2 , baseline calculation module  210  is shown in greater detail, according to an exemplary embodiment. Baseline calculation module  210  includes data clean-up module  212 . Data clean-up module  212  generally receives data from the BMS computer system of the building and pre-filters the data for data synchronization module  214  and the other modules of baseline calculation module  210 . Data clean-up module  212  includes outlier analysis module  256 , data formatting module  258 , and sorting module  260  for pre-filtering the data. Data clean-up module  212  uses sub-modules  256 - 260  to discard or format bad data by normalizing any formatting inconsistencies with the data, removing statistical outliers, or otherwise preparing the data for further processing. Data formatting module  258  is configured to ensure that like data is in the same correct format (e.g., all time-based variables are in the same terms of hours, days, minutes, etc.). Sorting module  260  is configured to sort data for further analysis (e.g., place in chronological order, etc.). 
     Outlier analysis module  256  is configured to test data points and determine if a data point is reliable. For example, if a data point is more than a threshold (e.g., three standard deviations, four standard deviations, or another set value) away from the an expected value (e.g., the mean) of all of the data points, the data point may be determined as unreliable and discarded. Outlier analysis module  256  may further calculate the expected value of the data points that each data point is to be tested against. Outlier analysis module  256  may be configured to replace the discarded data points in the data set with a NaN or another flag such that the new value will be skipped in further data analysis. 
     According to another exemplary embodiment, outlier analysis module  256  can be configured to conduct a cluster analysis. The cluster analysis may be used to help identify and remove unreliable data points. For example, a cluster analysis may identify or group operating states of equipment (e.g., identifying the group of equipment that is off). A cluster analysis can return clusters and centroid values for the grouped or identified equipment or states. The centroid values can be associated with data that is desirable to keep rather than discard. Cluster analyses can be used to further automate the data clean-up process because little to no configuration is required relative to thresholding. 
     Data clean-up module  212  may further include any other pre-filtering tasks for sorting and formatting the data for use by baseline calculation module  210 . For example, data clean-up module  212  may include an integrator or averager which may be configured to smooth noisy data (e.g., a varying number of occupants in a building area). The integrator or averager may be used to smooth data over a desired interval (e.g., a 15 minute average, hourly average, etc.). 
     Baseline calculation module  210  includes data synchronization module  214 . Data synchronization module  214  is configured to select a possible set of variables estimated to be significant to energy usage in the building. Data synchronization module  214  selects the possible set of variables (e.g., a preliminary set of variables) that are provided to stepwise regression module  284  for selection of the actual set of variables to use to generate the baseline model. According to various exemplary embodiments, the selection of some or all of the set of variables to use for baseline model generation may occur in data synchronization module  214 , stepwise regression analysis  284 , or a combination of both. Data synchronization module  214  includes sub-modules for calculating predictor variables and selecting one or more of the predicted variables to include in the possible set of variables. Data synchronization module  214  further includes sub-modules for selecting observed (e.g., measured) data points for the set of variables. 
     According to one exemplary embodiment, data synchronization module  214  is configured to calculate degree days and energy days (e.g., a predictor variable associated with heating or cooling of a building) and determine which of these predictors should be used to yield a better baseline model. The outputs of data synchronization module  214  (e.g., inputs provided to regression analysis module  216 ) may include the measurements or predictor variables to use, a period of time associated with the measurements or predictor variables, and errors associated with the data included in the measurements or predictor variables. 
     Data synchronization module  214  includes enthalpy module  262 , balance point module  264 , model determination module  266 , regression period module  268 , integration module  270 , NaN module  272 , missing days module  274 , workdays module  276 , and observed variable selection module  278 . Enthalpy module  262  is configured to calculate an enthalpy given a temperature variable and a humidity variable. Enthalpy module  262  combines an outdoor temperature variable and an outside air humidity variable via a nonlinear transformation or another mathematical function into a single variable. The single variable may then be used by baseline calculation module  210  as a better predictor of a building&#39;s energy use than using both temperature and humidity values separately. 
     Balance point module  264  is configured to find an optimal balance point for a calculated variable (e.g., a variable based on an enthalpy value calculated in enthalpy module  262 , an outdoor air temperature variable, etc.). Balance point module  264  determines a base value for the variable for which the estimated variance of the regression errors is minimized. Model determination module  266  is configured to determine a type of baseline model to use for measuring and verifying energy savings. The determination may be made based on an optimal balance point generated by balance point module  264 . Modules  264 ,  266  are described in greater detail in  FIGS. 4A-4E . 
     Regression period module  268  is configured to determine periods of time that can be reliably used for model regression by baseline calculation module  210  and data synchronization module  214 . Regression period module  268  may identify period start dates and end dates associated with calculated and measured variables for the data synchronization. Regression period module  268  may determine the start date and end date corresponding with the variable with the longest time interval (e.g., the variable for which the most data is available). For example, regression period module  268  determines the period by finding the period of time which is covered by all variables, and providing the start date and end date of the intersection to data synchronization module  214 . Regression period module  268  is further configured to identify data within the periods that may be erroneous or cannot be properly synchronized. 
     Integration module  270  is configured to perform an integration over a variable structure from a given start and end time period (e.g., a time period from regression period module  268 ). According to an exemplary embodiment, integration module  270  uses a trapezoidal method of integration. Integration module  270  may receive an input from balance point module  264  or another module of data synchronization module  214  for performing an integration for a balance point determined by balance point module  264 . NaN module  272  is configured to identify NaN flags in a variable structure. NaN module  272  is further configured to replace the NaN flags in the variable structure via interpolation. NaN module  272  may receive an input from, for example, data clean-up module  212 , and may be configured to convert the outlier variables and NaNs determined in module  212  into usable data points via interpolation. 
     Missing days module  274  is configured to determine days for which is there is not enough data for proper integration performance. Missing days module  274  compares the amount of data for a variable for a given day (or other period of time) and compares the amount to a threshold (e.g., a fraction of a day) to make sure there is enough data to accurately calculate the integral. Workdays module  276  is configured to determine the number of work days in a given interval based on the start date and end date of the interval. For example, for a given start date and end date, workdays module  276  can determine weekend days and holidays that should not figure into the count of number of work days in a given interval. Modules  274 ,  276  may be used by data synchronization module  214  to, for example, identify the number of days within a time interval for which there exists sufficient data, identify days for which data should not be included in the calculation of the baseline model, etc. 
     Observed variable selection module  278  is configured to receive observed or measured data from the BMS and determine which observed data should be used for baseline model generation based on the selection of calculated data in modules  264 - 266 . For example, when balance point module  264  determines a calculated variable, observed variable selection module  278  is configured to determine if there is enough predictor variable data for the observed variable. According to an exemplary embodiment, the predictor variable data and observed variable data for a specific variable (e.g., temperature) may only be used when sufficient predictor variable data (e.g., degree days) for the observed variable data exists. For example, if the predictor variable data is available over a specified range (e.g., 20 days, 2 months, or any other length of time), then module  278  may determine there is enough predictor variable data such that the predictor variable data and observed variable data can be used for baseline model generation. Observed variable selection module  278  is described in greater detail in  FIG. 3A . 
     Baseline calculation module  210  further includes regression analysis module  216 . Regression analysis module  216  is configured to generate the baseline model via a PLSR method. Regression analysis module  216  includes baseline model generation module  280  for generating the baseline model and PLSR module  282  for receiving data from data synchronization module  214 , applying the data to a PLSR method for, and providing baseline model generation module  280  with the method output. 
     Baseline model generation module  280  is configured to generate the baseline model. Baseline model generation module  280  is configured to use PLSR module  282  to perform PLSR of the data and stepwise regression module  284  to determine the predictor variables for the baseline model and to eliminate insignificant variables. Module  280  is configured to provide, as an output, the baseline model along with calculating various statistics for further analysis of the baseline model (e.g., computing the number of independent observations of data in the data set used, computing the uncertainty of the model, etc.). 
     Regression analysis module  216  is further shown to include stepwise regression module  284 . Stepwise regression module  284  is configured to perform stepwise linear regression in order to eliminate statistically insignificant predictor variables from an initial set of variables selected by data synchronization module  214 . In other words, stepwise regression module  284  uses stepwise regression to add or remove predictor variables from a data set (e.g., the data set from data synchronization module  214 ) for further analysis. 
     A stepwise regression algorithm of module  284  is configured to add or remove predictor variables from a set for further analysis in a systematic way. At each step the algorithm conducts statistical hypothesis testing (e.g., by computing a probability of obtaining a test statistic, otherwise known as a p-value, of an F-statistic, which is used to describe the similarity between data values) to determine if the variable should be added or removed. For example, for a particular variable, if the variable would have a zero (or near zero) coefficient if it were in the baseline model, then the variable is removed from consideration for the baseline model. According to various alternative embodiments, other approaches to stepwise regression are used (e.g., factorial designs, principal component analysis, etc.). Referring also to  FIG. 10 , instructions to add or remove variables from future consideration based on the analysis of module  216  may be provided to, for example, input calibration module  236  for affecting the queries run by data retriever  208 . 
     PLSR module  282  is configured to receive a subset of the variables from data synchronization module  214  which has been selected by stepwise regression module  284 , and to compute a partial least squares regression of the variables in order to generate a baseline model. According to various alternative embodiments, other methods (e.g., a principal component regression (PCR), ridge regression (RR), ordinary least squares regression (OLSR)) are also or alternatively used in the baseline model calculation instead of a PLSR method. 
     Baseline models calculated using historical data generally include four possible sources of error: modeling errors, sampling errors, measurement errors, and errors relating to multiple distributions in the data set. Sampling errors occur when the number of data samples used is too small or otherwise biased. Measurement errors occur when there is sensor or equipment inaccuracy, due to physics, poor calibration, a lack of precision, etc. Modeling errors (e.g., errors associated with the data set) occur due to inaccuracies and inadequacies of the algorithm used to generate the model. Errors relating to multiple distributions in the data set occur as more data is obtained over time. For example, over a one to three year period, data may be collected for the period and older data may become obsolete as conditions change. The older data may negatively impact the prediction capabilities of the current baseline model. 
     Conventional baseline energy calculations use ordinary least squares regression (OLS). For example, ASHRAE Guideline 14-2002 titled “Measurement of Energy Demand Savings” and “The International Performance Measurement and Verification Protocol” (IPMVP) teach that OLS should be used for baseline energy calculations. For OLS:
 
 y= 1*β 0   +Xβ   OLS +ε
 
where y is a vector of the response variables, X is a matrix consisting of n observations of the predictor variables, β 0  is an unknown constant, β OLS  is an unknown vector of OLS regression coefficients, and is a vector of independent normally distributed errors with zero mean and variance σ 2 . The regression coefficients are determined by solving the following equation:
 
β OLS =( X   T   X ) −1   X   T   y.  
 
     PLSR may outperform OLS in a building environment where the inputs to an energy consumption can be many, highly correlated, or collinear. For example, OLS can be numerically unstable in such an environment resulting in large coefficient variances. This occurs when X T X, which is needed for calculating OLS regression coefficients, becomes ill-conditioned in environments where the inputs are many and highly correlated or collinear. In alternative embodiments, PCR or RR are used instead of or in addition to PLSR to generate a baseline model. In the preferred embodiment PLSR was chosen due to its amenability to automation, its feature of providing lower mean square error (MSE) values with fewer components than methods such as PCR, its feature of resolving multicollinearity problems attributed to methods such as OLS, and due to its feature of using variance in both predictor and response variables to construct model components. 
     Baseline calculation module  210  is further shown to include cross-validation module  218 . Cross-validation module  218  is configured to validate the baseline model generated by regression analysis module  216  (e.g., there is no overfitting of the model, the order and number of variables in the model is correct, etc.) by applying data for a test period of time (in the past) to the model and determining whether the model provides a good estimate of energy usage. Cross-validation of the baseline model is used to verify that the model will fit or adequately describe varying data sets from the building. According to one exemplary embodiment, cross-validation module  218  may use a K-fold cross-validation method. The K-fold cross validation method is configured to randomly partition the historical data provided to baseline calculation module  210  into K number of subsamples for testing against the baseline model. In other embodiments, a repeated random sub-sampling process (RRSS), a leave-one-out (LOO) process, a combination thereof, or another suitable cross-validation routine may be used by cross-validation module  218 . 
     Referring now to  FIG. 3A , a flow chart of a process  290  for determining observed or measured variables to use in generation of a baseline model is shown, according to an exemplary embodiment. Process  290  is configured to select observed variables based on predictor variables generated by the data synchronization module of the baseline calculation module. Process  290  includes receiving data (step  291 ). Process  290  further includes determining the largest period of time for which there is data for predictor variables and observed variables (step  292 ). The period of time determined in step  292  may represent a period of time for which there will be enough predictor variable data for the corresponding data received in step  291 . Step  292  may include, for example, removing insufficient data points and determining the longest period for which there is enough data. For example, if there is too little data for one day, it may be determined that a predictor variable for that day may not be generated and therefore the day may not be used in ultimately determining a baseline model. 
     Process  290  includes initializing the observed variable (step  293 ). Initializing the observed variable includes determining a start and end point for the observed variable data, determining the type of data and the units of the data, and any other initialization step. Step  293  is used to format the received data from step  291  such that the observed data is in the same format as the predictor variable data. 
     Process  290  includes determining if enough predictor variable data exists (step  294 ). For example, if there is enough predictor variables (e.g., energy days) for a set period of time (e.g., 20 days), then process  290  determines that the predictor variables and its associated observed variable (e.g., enthalpy) may be used for baseline model generation. 
     Referring now to  FIG. 3B , a flow chart of a process  300  for determining calculated variables to use in generation of a baseline model is shown, according to an exemplary embodiment. Selecting some calculated variables for inclusion in a regression analysis used to generate a baseline model may provide better results than selecting some other calculated variables for inclusion, depending on the particulars of the building and its environment. In other words, proper selection of calculated variables can improve a resultant baseline model&#39;s ability to estimate or predict a building&#39;s energy use. Improvements to energy use prediction or estimation capabilities can improve the performance of algorithms that rely on the baseline model. For example, an improved baseline model can improve the performance of demand response algorithms, algorithms for detecting abnormal energy usage, and algorithms for verifying the savings of an energy conservation measure (e.g., M&amp;V calculations, etc.). 
     Process  300  provides a general process for selecting calculated variables to use in generation of a baseline model.  FIGS. 4A-4E  provide a more detailed view of process  300 . The output of process  300  (and of the processes shown in  FIGS. 4A-4E ) is the selection of calculated variables to use to generate the baseline model. Particularly, in an exemplary embodiment, process  300  selects between cooling energy days, heating energy days, cooling degree days, and heating degree days. The selection relies on a calculation of balance points (e.g., optimal base temperatures or enthalpies) of a building and using the calculations to calculate the potential variables (e.g., the energy days and degree days) for selection into the set of variables used to generate the baseline model. 
     The calculation and selection of calculated variables (for inclusion into the baseline model generation) is based in part on calculated balance points and may be accomplished in different ways according to different exemplary embodiments. According to one embodiment, a nonlinear least squares method (e.g., a Levenburg-Marquardt method) may be used to find the best calculated variables. Such a method, for example, may use daily temperature and energy meter readings to calculate balance points. A nonlinear least squares method may then be applied to the balance points to generate and select the appropriate calculated variables. 
     According to another embodiment, an optimization scheme may be used to determine the best balance point or points. The optimization scheme may include an exhaustive search of the balance points, a gradient descent algorithm applied to the balance points to find a local minimum of the balance points, a generalized reduced gradient method to determine the best balance point, and a cost function that is representative of the goodness of fit of the best balance point. The cost function may be an estimated variance of the model errors obtained from an iteratively reweighted least squares regression method, according to an exemplary embodiment. The iteratively reweighted least squares regression method is configured to be more robust to the possibility of outliers in the set of balance points generated and can therefore provide more accurate selections of calculated variables. 
     The optimization scheme algorithm may also use statistics (e.g., a t-statistic representative of how extreme the estimated variance is) to determine if building energy use is a function of, for example, heating or cooling. The statistics may be used to determine which balance points to calculate as necessary (e.g., calculating balance points relating to heating if statistics determine that building energy use is based on heating the building. 
     Referring to  FIG. 3B  and  FIGS. 4A-4E , an optimization scheme is described which uses a golden section search rule to calculate energy days and degree days and to determine which of the calculated variables to use based on a statistics to determine the type of energy use in the building. 
     Process  300  includes receiving data such as temperature data, humidity data, utility meter data, etc. (step  302 ). Process  300  further includes using the received data to calculate possible balance points (step  304 ). For example, step  304  may include using the received temperature data and humidity data to calculate an enthalpy. As another example, step  304  may include determining an optimal base temperature using the received temperature data. Process  300  further includes steps  306 - 318  for determining a calculated variable to use for baseline model generation based on enthalpy and temperature data calculated in step  304 ; according to various exemplary embodiments, steps  306 - 318  of process  300  may be used to determine calculated variables based on other types of balance points. 
     Process  300  includes steps  306 - 310  for determining optimal predictor variables based on the enthalpy calculated in step  304 . Process  300  includes determining a type of baseline model for cooling or heating energy days using the enthalpy (step  306 ). Process  300  further includes finding an optimal enthalpy balance point or points and minimum error variance for the resultant cooling and/or heating energy days (step  308 ). The optimal enthalpy balance point relates to, for example, a preferred base enthalpy of the building, and the minimum error variance relates to the variance of the model errors at the optimal balance point (determined using IRLS). Process  300  further includes determining if the optimal predictors determined in step  308  are significant (step  310 ). 
     Process  300  includes steps  312 - 316  for determining optimal predictor variables based on a temperature (e.g., temperature data received in step  302  by baseline calculation module  210 ). Process  300  includes determining a type of baseline model for cooling or heating degree days using the temperature (step  312 ). Process  300  further includes finding an optimal temperature balance point and minimum error variance for the cooling and/or heating degree days (step  314 ). Process  300  also includes determining if the optimal predictors determined in step  314  are significant (step  316 ). Using the results of steps  306 - 316 , process  300  determines which of energy days and degree days yields a better (e.g., more accurate) baseline model (step  318 ) when used by baseline calculation module  210 . 
     Referring now to  FIGS. 4A-4E , a detailed flow chart of process  300  of  FIG. 3B  is shown, according to an exemplary embodiment. Process  400  of  FIGS. 4A-4E  is shown using enthalpy and temperature to determine the balance points. The balance points are used to calculate the optimal degree or energy days predictor variable and in determining which calculated variables to use for baseline model generation. According to other embodiments, other methods may be used to determine the balance points. Referring more specifically to process  400  shown in  FIG. 4A , process  400  may calculate an enthalpy using temperature data input  402  and humidity data  404  (step  408 ). According to an exemplary embodiment, enthalpy may be calculated using a psychometric calculation. Process  400  includes receiving meter data  406  and enthalpy data and averages the data over all periods (step  410 ) for use in the rest of process  400 . 
     Process  400  includes determining possible baseline model types (i.e., whether both the heating and cooling balance points are needed to describe energy use in the building) based on the calculated enthalpy (step  412 ). For example, step  412  includes the method of determining a predictor variable associated with minimum energy use and then sorting all of the calculated variables (e.g., the variables determined in steps  408 - 410 ) and finding where the minimum energy predictor variable ranks compared to the other predictor variables. 
     Process  400  includes determining if using cooling base enthalpy in the baseline model calculation is feasible (step  414 ). If the predictor variable associated with the minimum energy found in step  412  is close enough to the maximum calculated variable, then it may be determined that a cooling base does not exist because cooling does not significantly impact energy consumption in the building or it cannot be found due to lack of data. If using the cooling base enthalpy is not feasible, the cooling base enthalpy is set to NaN and the minimum sigma is set to infinity (step  428 ) such that both values will be “ignored” later by process  400 . 
     If using a cooling base enthalpy is feasible, a range of feasible cooling base enthalpies is set (step  416 ). The range may vary from the maximum average monthly enthalpy to ten units less than the predictor variable associated with the minimum energy use. 
     Process  400  includes finding the base temperature of the predictor variable (e.g., via balance point module  264 ) by finding the base enthalpy for which the estimated variance of the regression errors is minimized. According to one exemplary embodiment, the minimization may be performed using the golden section search rule. Process  400  includes initializing the golden section search (step  418 ) and iterating the golden section search (step  420 ) until a desired base tolerance has been reached (step  422 ). The base tolerance may be predetermined via a logarithmic function of the size of the range, according to an exemplary embodiment. The golden section search of steps  420 - 422  provides an optimal balance point. The optimal balance point is then used to calculate a measure of variability and determine the t-statistic for the predictor variable. 
     When a desired base tolerance has been reached for the golden section search (step  422 ), process  400  may determine whether the t-statistic is significant (step  424 ). If the t-statistic is not significant, the minimum sigma representative of the t-statistic is set to infinity (step  428 ). If the t-statistic is significant, it is used in a later step of process  400  to determine the best predictor variable to use for the baseline model. 
     Referring now to  FIG. 5 , the objective function used in the golden section search of  FIGS. 4A-4E  is shown in greater detail, according to an exemplary embodiment. Process  500  is configured to calculate the objective function for use in the golden section search. Process  500  includes receiving data from, for example, step  416  of process  400  relating to a range of enthalpies or temperatures (or other measurements) that may be used for the baseline model. Process  500  includes, for all periods, finding an average predictor variable for each given balance point (step  502 ). For example, the following integral may be used to find the predictor variable: 
               1   T     ⁢       ∫   periodstart   periodend     ⁢       max   ⁡     (     0   ,       X   ⁡     (   t   )       -   b       )       ⁢   dt             
while the following integral may be used to determine the average response variable:
 
               1   T     ⁢       ∫   periodstart   periodend     ⁢       Y   ⁡     (   t   )       ⁢   dt             
where b is the balance point and T is the length of the period.
 
     After obtaining the predictor variable, process  500  includes performing an iteratively reweighted least squares method (IRLS) (step  504 ). IRLS is used because it is more robust to outliers than standard OLS methods. Process  500  includes using the results of step  504  to obtain an estimate of the error variance (step  506 ) which is used by process  400  to determine the predictor variable with the best fit for generating a baseline model. 
     Referring back to  FIGS. 4A-4B , process  400  further includes repeating the steps of steps  414 - 428  for the heating base enthalpy instead of the cooling base enthalpy. Referring now to  FIG. 4B , process  400  includes determining if heating base enthalpy is feasible (step  430 ), setting a range of feasible heating base enthalpies (step  432 ), initializing a golden section search (step  434 ), iterating the golden section search (step  436 ) until a desired base tolerance is reached (step  438 ), and determining if the t-statistic is significant (step  440 ). If the t-statistic is not significant, the minimum sigma is set to infinity (step  444 ), and if otherwise, the t-statistic will be used later in process  400 . 
     Process  400  further includes repeating the steps shown in  FIGS. 4A-B , only for the temperature instead of the enthalpy. Referring now to  FIG. 4C , process  400  includes determining possible model types based on the temperature data (step  446 ). Process  400  further includes determining if cooling base temperature is feasible (step  448 ), setting a range of feasible cooling base temperatures (step  450 ), initializing a golden section search (step  452 ), iterating the golden section search (step  454 ) until a desired base tolerance is reached (step  456 ), and determining if the t-statistic is significant (step  458 ). Referring now to  FIG. 4D , process  400  includes determining if heating base temperature is feasible (step  464 ), setting a range of feasible heating base temperatures (step  466 ), initializing a golden section search (step  468 ), iterating the golden section search (step  470 ) until a desired base tolerance is reached (step  472 ), and determining if the t-statistic is significant (step  474 ). If the t-statistic is insignificant for either, the cooling or heating base temperature is set to NaN and the minimum sigma for the cooling or heating base temperature is set to infinity (steps  462 ,  478  respectively). 
     Process  400  is then configured to recommend a predictor variable based on the base temperatures and minimum sigmas determined in the process. Process  400  includes recommending a default cooling degree day calculation (step  484 ) as a predictor variable if both the cooling base temperature and cooling base enthalpy were both set to NaN in process  400  (step  480 ). Process  400  may also recommend cooling degree days as a predictor variable if the minimum sigma for cooling energy days is better (e.g., lower) than the minimum sigma for cooling degree days (step  482 ). Otherwise, process  400  recommends using cooling energy days (step  486 ). 
     Process  400  may repeat steps  488 - 494  for heating degree days and heating energy days. Process  400  includes recommending a default heating degree day calculation (step  492 ) as a predictor variable if both the heating base temperature and heating base enthalpy were both set to NaN in process  400  (step  488 ). Process  400  may also recommend heating degree days as a predictor variable if the minimum sigma for heating energy days is better than the minimum sigma for heating degree days (step  490 ). Otherwise, process  400  recommends using heating energy days (step  494 ). 
     Referring now to  FIG. 6 , a flow chart of a process  600  of calculating enthalpy is shown, according to an exemplary embodiment. Process  600  includes receiving temperature and humidity data (step  602 ). Step  602  may further include identifying and removing humidity data points that are NaN, converting temperature data points to the correct format, or any other pre-processing steps. 
     Process  600  further includes, for each temperature data point, finding a corresponding humidity data point (step  604 ). For example, for a given time stamp for a temperature data point, step  604  includes searching for a corresponding time stamp for a humidity data point. According to an exemplary embodiment, a humidity data point with a time stamp within 30 minutes (or another period of time) of the time stamp of the temperature data point may be chosen as a corresponding humidity data point. Step  604  may further include searching for the closest humidity data point time stamp corresponding with a temperature data point time stamp. If a corresponding humidity data point is not found for a temperature data point, an enthalpy for the time stamp of the temperature data point is not calculated. 
     Process  600  further includes, for each corresponding temperature and humidity data point, calculating the enthalpy for the corresponding time stamp (step  606 ). The enthalpy calculation may be made via a nonlinear transformation, according to an exemplary embodiment. The calculation includes: converting the temperature data into a Rankine measurement and calculating the partial pressure of saturation using the below equation: 
             pws   =     exp   ⁢     {         C   1     T     +     C   2     +       C   3     ⁢   T     +       C   4     ⁢     T   2       +       C   5     ⁢     T   3       +       C   6     ⁢     T   4       +       C   7     ⁢     ln   ⁡     (   T   )           }             
where C1 through C7 are coefficients and T is the temperature data. The coefficients may be, for example, based on ASHRAE fundamentals. The enthalpy calculation further includes: calculating the partial pressure of water using the partial pressure of saturation:
 
             pw   =     H     100   *   pws             
where H is the relative humidity data. The enthalpy calculation further includes calculating the humidity ratio:
 
             W   =       0.621945   *   pw       p   -   pw             
where W is in terms of pounds water per pound of dry air. The enthalpy calculation further includes the final step of calculating the enthalpy in BTUs per pound dry air:
 
Enthalpy=0.24* T+W *(1061+0.444* T )
 
Once the enthalpy is calculated, the enthalpy is used rather than temperature data or humidity data in regression analysis to generate the baseline model (step  608 ).
 
     Detecting Changes in Energy Use in the Building to Support Measurement and Verification Systems and Methods 
     As described above, the predictor variables for use for the above-described baseline model may generally be calculated as described above by the equation: 
             x   =       ∫     t   1       t   2       ⁢       max   ⁡     (     0   ,       p   ⁡     (   t   )       -   b       )       ⁢   dt             
where x is the synchronized predictor variable, p(t) is the input predictor variable, t 1  and t 2  represent the start and end points for the period of time, and b is the balance point. The synchronized predictor variables are then used to develop the linear model:
 
 ŷ   k =β 0   t   k +β 1   x   1,k +β 2   x   2,k + . . . +β p   x   p,k  
 
where t k  is the time duration of the period and β q  is the q th  marginal energy cost per synchronized predictor variable x q,k . One or more of the systems or methods described above may be used to develop the linear model and to select appropriate predictor variables for the baseline model.
 
     Once the baseline model is developed, it is used, e.g., during a reporting period, to estimate the amount of energy that would have been used during the reporting period had the ECM not been performed. The avoided energy use can be described and calculated using:
 
energy avoidance=BE±RA±NRA−RPE
 
where BE is the baseline energy (e.g., average energy use during the baseline), RA are the routine adjustments (e.g., adjustments based on the predictor variables and typically included in the baseline model), NRA are the non-routine adjustments (e.g., adjustments based on (typically constant) factors related to energy usage that are not included in the baseline model), and RPE is the reporting period energy (e.g., the actual reported energy usage in the building). If there are no non-routine adjustments, the energy avoidance can be represented as:
 
energy avoidance= ŷ   k −RPE
 
where ŷ k  is the estimated energy use using the baseline model, i.e., the estimated energy use (BE) plus the routine adjustments based on the reporting period predictor variables (RA).
 
     If the need for non-routine adjustments is identified (e.g., if static factors have changed in a way such that energy use is affected), then non-routine adjustments must be performed by taking them into account in the energy avoidance calculation. Changes in static factors that could necessitate non-routine adjustments may include, for example, changes in the size of the used facility space, the hours of operation of the building, the number of employees using the building, and the number of computer servers in the building. Any factor not expected to normally change and having an impact on energy use that is not included in the predictor variables of the baseline model could be considered a static factor. A change in static factors may occur during the reporting period or during the baseline period. 
     Referring generally to the figures below, systems and methods are described for automatically detecting changes in static factors during a baseline period or reporting period. A method for automatically detecting changes in static factors may include developing a null hypothesis of a baseline model with constant coefficients (stationary). The method further includes determining whether a statistically significant change in the baseline model&#39;s parameters has occurred. If statistically significant changes in the baseline model&#39;s parameters are detected over time, the method can include outputting an indication (e.g., to a user interface, to a computerized module for recalculating the baseline model, etc.) that a static factor has changed. The method can include conducting multiple sequential hypothesis tests spread over a period of time. Known theoretical correlations between temporally adjacent test statistics can be used to provide false positive suppression (e.g., suppress an indication that a static factor has changed where the data causing such an indication was spurious due to weather or another transient condition) without the undesirable effect of greatly reducing the test power (as would occur if, for example, the Bonferroni correction was applied). 
     When a building is operating in a consistent manner (e.g., consistent energy usage) and the baseline model for the building includes all the independent predictor variables necessary to accurately estimate the energy usage, the coefficients of the baseline model should remain constant over time. Therefore, if two temporally consecutive windows of data from time intervals [t a ,t b ] and [t b ,t c ] are used, the difference in two baseline model coefficients should be near zero. The difference in model coefficients can be represented as:
 
Δβ={circumflex over (β)} 1 −{circumflex over (β)} 2  
 
where Δβ is the difference between the baseline model coefficients from window one and window two {circumflex over (β)} 1 , {circumflex over (β)} 2 , respectively. Because the baseline model coefficients have physical meaning (e.g., cost per cooling degree day), unexpected changes in coefficients over time can advantageously be linked to root causes (e.g., chiller fouling, decrease in setpoints, etc.).
 
     For the coefficients, there may be random variation in a coefficient the magnitude of which is based on, for example: the number of periods or data points in the time intervals, the variance of the errors of the baseline model, the number of predictor variables used in the model, and the values of the predictor variables during each of the two time intervals. Additionally, the values of the predictor variables during each time interval can have a significant effect of the variation of the coefficients. Thus, in the preferred embodiment, the statistical methods described below may be used to determine whether the difference in the coefficients is large enough to be considered statistically significant or whether the coefficient difference is due to the random variation described above, rather than a real change in a static factor affecting the building&#39;s energy use. 
     Referring now to  FIG. 7 , a flow chart of a process  700  for automatically detecting changes in static factors that occurred between two sequential time periods is shown, according to an exemplary embodiment. Process  700  includes developing a null hypothesis that the baseline model is accurate (e.g., the baseline model is a correct representation of building energy usage) and stationary (e.g., the coefficients are constant) (step  702 ). 
     Process  700  further includes calculating the baseline model including the baseline model coefficients for two different chronologically ordered data sets (step  704 ). The coefficients may include, for example, energy use per degree day or energy use per energy day. In other embodiments, the coefficients could include an aggregation, or other value that describes the coefficients for a given periodic baseline model. In step  705 , process  700  includes calculating at least one test statistic that is related to the difference in a coefficient or a set of coefficients from the two data sets. If two consecutive windows of data are used to build similar baseline models (i.e., the coefficients of the models are similar) and static factor changes have not occurred during the time period of the windows, then the test statistic should be small (i.e., within the expected amount of random variation). 
     Process  700  further includes determining a probability density function (PDF) of the test statistic (step  706 ). Using the PDF of the test statistic, process  700  determines a critical value such that the probability of the test statistic being greater than the critical value is the same as the specified probability of falsely rejecting the null hypothesis of step  702  (step  708 ). A user of the building management system may select an acceptable level for the probability of falsely rejecting the null hypothesis, according to an exemplary embodiment. In some embodiments, the probability may be built into contract standards and differ from building to building, building type to building type, or otherwise vary. 
     Process  700  further includes rejecting or accepting the null hypothesis based on whether the test statistic is greater than the critical value found in step  708  (step  710 ). If the test statistic is greater than the critical value, then the null hypothesis of a constant baseline model is rejected. When the null hypothesis is rejected, the system can determine that a static factor has changed and a non-routine adjustment may be necessary. The user may be alerted to the necessity of a non-routine adjustment, additional calculations may be performed, or other various actions related to the non-routine adjustment may be performed by the building management system. 
     Referring further to process  700 , least squares estimation may be used to find the coefficients of the baseline models. The least squares estimation problem can be stated as the following: given a linear model
 
 Y=Xβ+ε, ε˜N (0,σ 2   I )
 
find the vector {circumflex over (β)} that minimizes the sum of squared error RSS:
 
RSS=∥ Y−X{circumflex over (β)}∥   2 .
 
     In the above equations, Y is a vector that contains the individual n observations of the dependent variable and X is a n by p+1 matrix that contains a column of ones and the p predictor variables at which the observation of the dependent variable was made. ε is a normally distributed random vector with zero mean and uncorrelated elements. According to various exemplary embodiments, other methods than using RSS may be used (e.g., weighted linear regression, regression through the origin, etc.) The optimal value of {circumflex over (β)} based on the least squares estimation has the solution stated above with respect to  FIG. 2 :
 
{circumflex over (β)}=( X   T   X ) −1   X   T   Y  
 
and {circumflex over (β)} is a normal random vector distributed as:
 
{circumflex over (β)}˜ N (β,σ 2 ( X   T   X ) −1 ).
 
The resulting sum of squared error divided by sigma squared is a chi-square distribution:
 
     
       
         
           
             
               
                 RSS 
                 
                   σ 
                   2 
                 
               
               ~ 
               
                 𝒳 
                 
                   n 
                   - 
                   
                     ( 
                     
                       p 
                       + 
                       1 
                     
                     ) 
                   
                 
                 2 
               
             
             . 
           
         
       
     
     The difference in the coefficients is distributed as
 
Δβ={circumflex over (β)} 1 −{circumflex over (β)} 2   ˜N (0,σ 2 [( X   1   T   X   1 ) −1 +( X   2   T   X   2 ) −1 ]).
 
The quadratic form of a normally distributed random vector where the symmetric matrix defining the quadratic form is given by the inverse of the covariance matrix of the normal random vector is itself a chi-square distributed random variable with degrees of freedom equal to the length of Δβ:
 
                         Δβ   T     ⁡     [         (       X   1   T     ⁢     X   1       )       -   1       +       (       X   2   T     ⁢     X   2       )       -   1         ]         -   1       ⁢   Δβ       σ   2       ~     𝒳     p   +   1     2       .         
Additionally, the sum of two independent chi-square distributions is itself a chi-square distribution with degrees of freedom equal to the sum of the degrees of freedom of the two original chi-square distributions. Thus, the sum of the two sum of squared errors divided by the original variance is chi-square distributed, as:
 
                     RSS   1     +     RSS   2         σ   2       ~     𝒳       n   1     +     n   2     -     2   ⁢     (     p   +   1     )         2       .         
n 1  and n 2  are the number of data points used to estimate the model coefficients {circumflex over (β)} 1 , {circumflex over (β)} 2 . Finally, the ratio of two chi-square distributions divided by their respective degrees of freedom is an F-distributed random variable:
 
               F   Δβ     =       (             Δβ   T     ⁡     [       (       X   1   T     ⁢     X   1       )     +       (       X   2   T     ⁢     X   2       )       -   1         ]         -   1       ⁢   Δβ         RSS   1     +     RSS   2         )     ⁢         (         n   1     +     n   2     -     2   ⁢     (     p   +   1     )           p   +   1       )     ~     F       p   +   1     ,       n   1     +     n   2     -     2   ⁢     (     p   +   1     )               .             
F Δβ  is defined as the test statistic. As Δβ moves away from the origin, F Δβ  increases. Further, the maximum increase occurs in the direction of the least variance of the model coefficients and is scaled by the sum of squared errors. Thus, F Δβ  is based on changes in model coefficients which can easily be related back to a root cause and it takes into account the random variation of the changes of the model coefficients even when the model is stationary. The F Δβ  statistic may further be converted into a standard normal variable Z Δβ  by the proper transformation function.
 
     The F Δβ  or statistic can be used as the test statistic in step  710  of process  700  (if Z Δβ  is used then the critical value will be calculated with the inverse distribution function of the standard normal distribution). The user of the building (or a contract or automated process) determines an acceptable level for α, the probability of rejecting the null hypothesis when it is in fact valid. An automated process uses α to determine the critical value for use in accepting or rejecting the null hypothesis. The null hypothesis is rejected if the F-statistic F Δβ  is greater than its critical value f crit  which may be calculated using F p+1,n     1     +n     2     −2(p+1)   −1 (1−α), where F −1  is the inverse of the cumulative F-distribution with the required degrees of freedom. In other words, the null hypothesis is rejected and a static factor can be determined to have changed when F Δβ  &gt;f crit . 
     Process  700  is used to determine whether a single test statistic, from a single test, including the model coefficients of two (time adjacent) building models, indicates a change in static factors in a building. Process  700  is more particularly used to determine whether a change has occurred at any point in time during the baseline period or reporting period. In an exemplary embodiment, process  700  is repeated multiple times over sets of data that are windowed using different methods based on whether the baseline period data (model building) or reporting period data (facility monitoring) is being tested for static factor changes (i.e., model coefficient consistency). Two representative different methods are shown in greater detail in  FIGS. 8A and 9A . In these cases, the critical value must be adjusted from that shown above to limit the false alarm rate of any of the multiple tests performed. 
     When attempting to detect static factor changes during the baseline or reporting period, several iterations of model calculations and associated coefficient testing can be conducted. For example, over the baseline data, the test may be run n−2(p+1)+1 times, where n is the number of data points and p is the number of predictor variables in each model. If that many tests are run, the probability of finding at least one false alarm will be much greater than a (in fact, as the number of tests increase, the probability approaches 1). Therefore, for the baseline period, some embodiments limit the probability that any of the several tests falsely reject the null hypothesis to avoid false alarms (i.e., false indications that a static factor of a baseline model has changed), since several tests will be performed. For the reporting period, limiting the probability of a false alarm in a set time duration to a predefined value can help avoid falsely rejecting the null hypothesis. For example, the user or electronics (e.g., based on minimum performance or contract specification) may specify that the probability of a false rejecting the null hypothesis on any test performed in a given year should be less than 15%. 
     Referring now to  FIG. 8A , a flow chart of a process  800  for detecting changes in static factors of a baseline model during a baseline period is shown, according to an exemplary embodiment. To detect changes, the Z Δβ  statistic is calculated at each possible value for the dividing time between the two data windows, and all baseline data (e.g., all data from the beginning of the baseline period to the beginning of the installation period) is used. That is, the statistic is calculated for each possible windowing (for which the data remains in chronological order and there is at least one degree of freedom in each window) in the baseline period. 
     Process  800  includes temporally windowing the data into two windows (step  802 ). During the baseline period, the data received and used in the baseline model is received by process  800  and divided into two different time periods (i.e., windows). The data is placed into the data windows in chronological order (e.g., into the first window until the first window fills, then the second window). Referring also to process  850  of  FIG. 8B , a process for determining the size and time periods for the data windows may be used. Baseline model coefficients may be calculated for the two different data windows, and the difference in coefficients Δβ is calculated (step  804 ). Step  804  may further include storing the difference (e.g., for consideration in step  820  of process  800 ). Process  800  further includes calculating the test statistic F Δβ  based on the change in model coefficients (step  806 ). Process  800  further includes converting the test statistic F Δβ  into a standard normal random variable Z Δβ  (step  808 ). The conversion is done by passing the test statistic through its cumulative probability distribution and then through the inverse of the standard normal distribution. The resultant standard normal random variable Z Δβ  is given by:
 
 Z   Δβ   =N   0,1   −1   {F   p+1,n     1     +n     2     −2(p+1) ( F   Δβ,k )}.
 
     Process  800  further includes determining whether to slide the data window (step  810 ). After calculating a difference in step  804  for two given time periods, the data windows for the two given time periods may be adjusted. For example, for two time periods defined by [t a , t b ] and [t b ,t c ], step  810  of sliding the data window is executed by increasing t b , increasing the first data window by one data point and decreasing the second data window by one data point. Sliding the data window may include moving data points from one time period to another based on a timestamp or other property of each individual data point. By sliding the data window (e.g., by moving the end point of the first data window and the start point of the second data window), new coefficients for the data windows may be calculated, and another difference in coefficients and Z-statistic may be calculated. Process  800  repeats steps  802 - 808  until, for example, the second data window contains the fewest allowable data points (e.g., the second window has p+2 data points). For example, process  800  may repeat steps  802 - 808 , beginning from when the first data window contains the fewest number of possible data point until the second data window contains the fewest number of possible data points. The method of repeating steps  802 - 808  by shifting the data windows is shown in greater detail in  FIG. 8B . 
     Referring now to  FIG. 8B , the shifting of the data windows for calculating a new observed Z Δβ  statistic is shown in greater detail. The detailed process of  FIG. 8B  may be the steps that occur during steps  802 - 806  of  FIG. 8A . Process  850  includes receiving response variables Y (e.g., variables representative of building energy usage)  852  and building predictor variables X  854 . Process  850  includes setting an initial first window end of p+2, where p is the number of predictor variables (step  856 ). Such a choice means that the initial first data window only includes p+2 data points (enough for one degree of freedom in the probability distribution calculations). The first window end marks the last data point that is in the first window; thus, the next data point is the first datum in the second window. If the observed maximum Z Δβ  statistic indicates that a static factor has changed, this point can be inferred to represent a point at or after which a static factor change occurred. For example, when coefficients are calculated using the two data windows and a difference in coefficients is large enough, the data point that is used as the endpoint or start point for the data windows may be determined to represent a time at or after which a potential change in static factors occurred. In alternative embodiments, the different windows could include gaps between the windows (e.g., data points “in between” the data points in the two windows) or some amount of overlap. 
     Process  850  further includes, while the first window end is less than or equal to n−p+2, increasing the first window end by one (e.g., sliding the first data window endpoint by one) (step  858 ). n is the number of data points in the entire baseline period being used. The data is divided into the two data windows based on the endpoints (step  860 ). Coefficients are calculated for each of the data windows (steps  862 ,  864 ). In step  862 , coefficients for the response variables Y for the two data windows are calculated, and in step  864 , coefficients for the predictor variables X for the two data windows are calculated. The difference in coefficients and the Δβ statistic is calculated (step  866 ), the F-statistic F Δβ  is calculated (step  868 ), and the standardized statistic Z Δβ  is calculated (step  870 ). Steps  866 - 870  corresponds with steps  804 - 808  of process  800 . Process  850  is then repeated until the second data window is eventually reduced to only p+2 data points. 
     Referring again to  FIG. 8A , after calculating the statistics, process  800  includes determining the maximum Z-statistic from steps  802 - 808  (step  812 ). Step  812  includes comparing the differences calculated in step  808  (e.g., or step  870 ) to find the maximum difference. Process  800  further includes calculating the P-value of the maximum Z-statistic (step  814 ). The P-value represents the probability of obtaining a maximum test statistic as extreme as the observed maximum Z-statistic determined in step  812  given that the null hypothesis (linear model with constant coefficients) is true. 
     The P-value may be calculated in various ways. In one exemplary embodiment, the P-value may be approximated using a Monte Carlo simulation. In one embodiment, Monte Carlo samples of the maximum Z Δβ  determined in step  812 . For example, steps  808 - 812  are performed on computer generated data that fits the null hypothesis, samples which have a maximum value greater than the data points used in process  800  may be determined, and the fraction or ratio of such samples may be used to calculate the P-value. In another exemplary embodiment, the P-value may be approximated using a Monte Carlo simulation of a multivariate normal vector. One exemplary use of a Monte Carlo simulation to calculate a P-value (e.g., step  814 ) is shown in greater detail in  FIG. 13 . 
     Process  800  further includes displaying the P-value to a user and displaying alarms or warnings based on the P-value if necessary (step  816 ). The P-value is representative of the likelihood of a change in the static factors of the building (e.g., a P-value near zero indicated that a change may have occurred). Step  816  includes using the P-value to generate an alarm or warning representative of the likelihood of the change. Step  816  may also or alternatively include providing the P-value to another system component or automated process (e.g., recalculating the baseline model with knowledge of the changed static factor likelihood). The energy cost or savings is determined if a change occurred (step  818 ). An adjustment to account for the change in static factors is also determined (step  820 ). Steps  818 - 820  may include estimations of energy cost or savings (e.g., based on the changed static factors and/or updated model) and adjustments based on the P-value. 
     Referring now to  FIG. 9A , a flow chart of a process  900  for detecting changes in static factors of a baseline model during a reporting period (as opposed to during a baseline period as described above with reference to  FIGS. 8A-B ) is shown, according to an exemplary embodiment. The reporting period may last as long as twenty years or more and several changes in static factors may occur over the time frame. Process  900  may be used once enough data is collected by the building management system, and after enough data is collected, old data may be discarded from process  900  as new data continues to be collected. 
     Process  900  includes calculating a critical value for the test statistic at a given alpha value (step  902 ). The alpha value can be supplied by the user or obtained from the standards of the contract and may vary from building to building or building type to building type. According to one embodiment, step  902  includes using a Monte Carlo cumulative probability distribution inversion to calculate the critical value. An exemplary method and use of such an inversion to find the critical value is shown in greater detail in  FIG. 12 . The samples used in the inversion may be performed using the actual data points (formed by simulating a linear model and calculating Z Δβ ) or may be simulated from a multivariable normal distribution. 
     Process  900  further includes aligning the two data windows with the last change (step  904 ). Step  904  includes receiving new data and sliding the data in the two windows. Data for the building during the reporting period is collected until two data windows are filled (step  906 ). In an exemplary embodiment, a window may not be considered full until a minimum amount of time has elapsed (e.g., one month, one year, etc.). The data in each window is then shifted forward (step  908 ) (e.g., as described in process  950  of  FIG. 9B ) and a normalized test statistic for the two data windows is calculated (step  910 ). Step  910  may include, for example, calculating coefficients for the two data windows and determining the difference in the coefficients as described in process  700  of  FIG. 7 . 
     Process  900  further includes determining whether the test statistic calculated in step  910  is greater than the critical value determined in step  902  (step  912 ). If the test statistic is greater than the critical value, the user may be alerted to a change in static factors or an indication of the change in static factors may be communicated (step  914 ) to another module of the building management system (e.g., a reporting module, a logging module, a module for automatically initiating recalculation of the energy model, etc.) and an energy cost or savings is determined (step  916 ). An adjustment to account for the change is determined as well (step  918 ). After the change or adjustment has completed, building response and predictor variables may be collected again until the initial windows are filled again. In other words, process  900  may repeat accordingly to continuously evaluate static factor changes in the utilized baseline model. 
     If the difference in step  910  is less than the critical value, process  900  repeats steps  908 - 910  by shifting the two data windows. The method of shifting the two data windows is shown in greater detail in  FIG. 9B . 
     Referring now to  FIG. 9B , the shifting of the data windows for calculating an new observed standardized statistic of step  908  of process  900  is shown in greater detail. Process  950  is configured to keep both data windows used at a fixed length, and moves the data windows forward in time instead of changing the size of either window. Therefore, when a new reporting period data point is obtained and used, the data point will become the last data point in the second data window, and the oldest data point in the first data window winds up being eliminated from the analysis. 
     Process  950  includes receiving a first window size, a second window size n 2 , and endpoints n 1 , n 2  (step  952 ). Process  950  includes using the endpoints for the two data windows (step  954 ). The first window size and second window size are the sizes of the data windows that are preferred for process  900 . The two sizes of the data windows will remain fixed throughout process  900 , and process  950  includes changing the endpoints of the data windows instead of the size of the data windows. The sizes of the data windows may be predetermined by the building management system or a user thereof. The endpoints n 1  and n 2  are used to determine the initial endpoints of the two windows (e.g., the first window includes all data points from a cutoff point to point n 1 , and the second window includes all data points from n 1 +1 to n 2 . 
     While the second window end is less than the reporting end (e.g., the end data point for the reporting period), the first window end n 1  and second window end n 2  are increased by one (step  956 ). The data received by process  950  is then divided into two windows with start points and endpoints defined by the first window end and second window end (step  958 ). Since the size of the two data windows are limited (e.g., based on the first window size and second window size specified in step  952 ), the “oldest” data point that used to be in the first data window in the previous iteration of process  950  no longer fits in the data window. Therefore, the data points in the first data window may range from a value z to the new n 1 . The second data window may continually be adjusted such that it includes data points from the new n 1 +1 to the new n 2 . 
     Coefficients are calculated for each of the data windows (steps  960 ,  962 ). In step  960 , the response variables are placed into two vectors corresponding to the two windows, and in step  962 , the predictor variables are placed into two matrices corresponding to the two windows. The difference in coefficients (and the Z Δβ  statistic) is calculated (steps  964 - 968 ). Process  950  is then repeated as new data is gathered throughout the reporting period until the second data window reaches the reporting period end. The Z Δβ  statistic is then used in the comparison of step  912  of process  900 . 
     Referring now to  FIG. 10 , a block diagram of baseline calculation module  210  is shown in greater detail, according to an exemplary embodiment. Baseline calculation module  210  is shown to include the data clean-up module, data synchronization module, regression analysis module, and cross-validation module as described in  FIG. 2 . Baseline calculation module  210  further includes data window module  1002 , coefficient calculation module  1004 , and test statistic module  1006 . Baseline calculation module  210  is configured to, using modules  1002 - 1006 , detect changes in static factors during a baseline period. Baseline calculation module  210  further includes predictor variables  1010  which may be provided to the various modules of the computer system (e.g., facility monitoring module  230 ). 
     Data window module  1002  is configured to perform the data window sliding and calculations as described in  FIGS. 8B and 9B . Data window module  1002  may receive data and either receive or determine initial data windows for the process of  FIG. 8B  to use. Data window module  1002  further slides the data windows (e.g., via processes  850 ,  950 ) and organizes the data in the two windows to facilitate the completion of the processes of  FIG. 8B . 
     Coefficient calculation module  1004  is configured to calculate the coefficients for the baseline model and for the detection of changes in static factors. Coefficient calculation module  1004  may be configured to, for example, conduct the calculations shown in steps  862  and  864  of  FIG. 8B  and steps  960  and  962  of  FIG. 9B . 
     Test statistic module  1006  is configured to receive coefficient calculation results from module  1004  and to calculate a test statistic Z Δβ  representative of the difference in coefficients. The calculation may include the calculation of the F-statistic F Δβ  which is then converted into the test statistic Z Δβ . Test statistic module  1006  may further be configured to calculate a PDF of the test statistic Z Δβ  and a PDF of the maximum of n test statistics Z Δβ . 
     Baseline calculation module  210  further includes energy savings module  1008 . Energy savings module  1008  is configured to determine the energy (or monetary) costs or savings associated with a change in static factors. 
     One method of estimating the cost or savings is to take data from the last detected change (or the end of the installation period) to the time the change was detected and the data available since the change occurred. Model coefficients can be calculated for both the pre- and post-change data. The estimated cost is then given by:
 
cost= x   T [β new −β old ]
 
with a standard error of:
 
 se (cost)={circumflex over (σ)}{1+ x   T [( X   old   T   X   old ) −1 +( X   new   T   X   new ) −1 ] x}   1/2 .
 
where β new ,β old  are the calculated coefficients, X old  is data from before the change, and X new  is data from after the change.
 
     Referring now to  FIG. 11 , a detailed block diagram of the facility monitoring module  230  of  FIG. 10  is shown, according to an exemplary embodiment. Facility monitoring module  230  is configured to detect changes in static factors during the reporting period. Module  230  receives baseline model information (e.g., predictor variables  1010 ) from baseline calculation module  210  and other building information such as false alarm probability data  1120  and other building data. False alarm probability data  1120  may be a user entered value that represents a desired probability of generating a false alarm of determining that a static factor has changed when it actually has not. For example, 5% may be a desired probability. False alarm probability data  1120  may further be received from system feedback module  234 , baseline agreement information  232 , or otherwise. 
     Facility monitoring module  230  includes critical value module  1102 . Critical value module  1102  is configured to calculate a critical value for the Z Δβ  statistic. Critical value module  1102  may be configured to perform a cumulative probability distribution inversion calculation as described in  FIG. 12 , according to an exemplary embodiment. 
     Facility monitoring module  230  further includes data window module  1104 , data collection module  1106 , and test statistic module  1108 . Data window module  1104  is configured to determine initial data window parameters for the processes of  FIGS. 8B and 9B  and may have the same functionality as data window module  1002  of  FIG. 10 . Data collection module  1106  collects building data and provides the data to data window module  1104  for filling the two data windows. Test statistic module  1108  is configured to calculate a test statistic using predictor variables  1010 . Test statistic module  1108  may calculate the test statistic may have the same functionality as test statistic module  1006  of  FIG. 10 . 
     Facility monitoring module  230  further includes energy savings module  1110  and reporting module  1112 . After calculating the test statistic and determining if a change in static factors has occurred, energy savings module  1110  is configured to calculate the energy savings or cost associated with the change, and reporting module  1112  is configured to generate a report about the change. Facility monitoring module  230  may use energy savings module  1110  and reporting module  1112  to generate a report to provide to system feedback module  234  for user output and analysis. Energy savings module  1110  may have the same functionality as energy savings module  1008  of  FIG. 10 , according to an exemplary embodiment. 
     Referring now to  FIG. 12 , process  1200  is used to find a critical value for a maximum of a sequence of standardized F-statistics (Z Δβ  statistics) where the sequence is determined using the process described in the present disclosure. Process  1200 , according to an exemplary embodiment, converts the problem of finding this critical value to a problem of finding a critical value for a maximum element of a multivariate normal vector with a known covariance matrix. In an exemplary embodiment, this problem is solved by performing a Monte Carlo simulation to find the inverse cumulative distribution function (CDF). Process  1200 , therefore, evaluates an inverse of a CDF to ultimately use the values to reject or accept a null hypothesis of an accurate and constant baseline model, according to an exemplary embodiment. 
     In multiple hypothesis tests, correlation of the test statistics may largely impact the conservativeness of typical methods for suppressing the family-wise probability of falsely rejecting the null hypothesis (Bonferroni&#39;s method, for example). For example, if multiple statistics two data sets are highly correlated, the statistics do not differ by a significant amount. Thus, direct applications of Bonferroni&#39;s method would be very conservative and greatly reduce the power of the test (probability of correctly identifying a change in the model coefficients). In the embodiments of the present disclosure, if static factors are not changing, the statistics calculated using the windowing method described previously should be highly correlated. Window data selection steps described above could be designed to maintain this high correlation during normal behavior. For example, during the reporting period, the last data point inserted into the second data window replaces the oldest data point in the first data window, meaning that only two data points have changed since the last calculation. Using the inverse CDF process of  FIG. 12 , values that are generated are suitable for use for determining the accuracy and consistency of a baseline model without having the issue of highly correlated data skewing the results. This is done by determining the autocorrelation of the test statistics, which is dependent on the method used to window the data and to slide the window forward, and using the correlation to directly calculate (or approximate) the inverse CDF of the maximum of this autocorrelated sequence. 
     Evaluation of the inverse CDF can be phrased as, given a value p (e.g., a desired probability), find a value x such that P(X&lt;x)=p, where X is a random variable, in the current disclosure the maximum of the sequence of statistics and x is the argument of the CDF, which in the current disclosure corresponds with the critical value of the null hypothesis (e.g., the value of step  708  of process  700  for the null hypothesis test). In context of the present disclosure, this means the inverse CDF is used to determine a critical value such that the probability that the maximum of the sequence of statistics is equal to the desired probability p typically equal to one minus the indicated probability of falsely rejecting the null hypothesis. 
     If several samples from drawn from the distribution of data points, a point estimate for the probability p is given by: 
     
       
         
           
             
               
                 P 
                 ^ 
               
               ⁡ 
               
                 ( 
                 
                   X 
                   &lt; 
                   x 
                 
                 ) 
               
             
             = 
             
               
                 p 
                 ^ 
               
               = 
               
                 
                   n 
                   
                     { 
                     
                       X 
                       &lt; 
                       x 
                     
                     } 
                   
                 
                 n 
               
             
           
         
       
     
     with an associated 1−α confidence interval for p given by the equation in step  1210  below. The confidence interval 1−α indicates a desired probability that the true value of p resides within the band {circumflex over (p)} plus or minus the tolerance. 
     Process  1200  includes determining a desired probability p (e.g., the p value of P&lt;=p) and a confidence interval for p (step  1202 ). The desired probability p and confidence interval may be chosen by a user of the building. The confidence interval should be determined such that probabilities with values on the upper and lower limits of the interval are accepted at the 1−α confidence level. Process  1200  is configured to return a value such that all probabilities within the 1−α confidence interval are included in the range defined by the upper and lower limits. This guarantees that the probability that the actual value of p for the returned value is between the upper and lower limits is greater than 1−α. 
     Process  1200  further includes determining the number of samples required to draw from the distribution in order to reach a desired tolerance (step  1204 ). The number of samples n may be found by using an iterative root finding technique where the objective is to find n such that: 
                 max   ⁡     (           n   ⁢     p   ^     ⁢       F       2   ⁢   n   ⁢     p   ^       ,     2   ⁡     [       n   ⁡     (     1   -     p   ^       )       +   1     ]           -   1       ⁡     (     a   /   2     )             n   ⁡     (     1   -     p   ^       )       +   1   +     n   ⁢     p   ^     ⁢       F       2   ⁢   n   ⁢     p   ^       ,     2   ⁡     [       n   ⁡     (     1   -     p   ^       )       +   1     ]           -   1       ⁡     (     a   /   2     )             -     low   ⁢           ⁢   limit       ,     
     ⁢       high   ⁢           ⁢   limit     -         (       n   ⁢     p   ^       +   1     )     ⁢       F       2   ⁡     [       n   ⁢     p   ^       +   1     ]       ,     2   ⁢     n   ⁡     (     1   -     p   ^       )             -   1       ⁡     (     1   -     a   /   2       )             n   ⁡     (     1   -     p   ^       )       +       (       n   ⁢     p   ^       +   1     )     ⁢       F       2   ⁡     [       n   ⁢     p   ^       +   1     ]       ,     2   ⁢     n   ⁡     (     1   -     p   ^       )             -   1       ⁡     (     1   -     a   /   2       )                 )       =   0     ,         
where {circumflex over (p)} is given the value of p and the low limit and high limit are the upper and lower limits of the 1−α confidence interval.
 
     Process  1200  further includes drawing n samples of the distribution at random (step  1206 ). The samples can be drawn by simulating a linear model and performing the process or in order to do an approximation from a multivariate normal distribution. Using the samples, a critical value x is found such that the total number of samples n drawn less than x is equal to np (e.g., the number of samples times the probability of each individual sample being less than x) and the total number of samples greater than x is equal to n(1−p) (e.g., the number of samples times the probability of each individual sample being greater than x) (step  1208 ). 
     Process  1200  further includes recalculating the 1−α confidence interval for p (step  1210 ). The equation used for the calculation may be the following: 
     
       
         
           
             
               P 
               ⁡ 
               
                 ( 
                 
                   
                     
                       n 
                       ⁢ 
                       
                         p 
                         ^ 
                       
                       ⁢ 
                       
                         
                           F 
                           
                             
                               2 
                               ⁢ 
                               n 
                               ⁢ 
                               
                                 p 
                                 ^ 
                               
                             
                             , 
                             
                               2 
                               ⁡ 
                               
                                 [ 
                                 
                                   
                                     n 
                                     ⁡ 
                                     
                                       ( 
                                       
                                         1 
                                         - 
                                         
                                           p 
                                           ^ 
                                         
                                       
                                       ) 
                                     
                                   
                                   + 
                                   1 
                                 
                                 ] 
                               
                             
                           
                           
                             - 
                             1 
                           
                         
                         ⁡ 
                         
                           ( 
                           
                             a 
                             / 
                             2 
                           
                           ) 
                         
                       
                     
                     
                       
                         n 
                         ⁡ 
                         
                           ( 
                           
                             1 
                             - 
                             
                               p 
                               ^ 
                             
                           
                           ) 
                         
                       
                       + 
                       1 
                       + 
                       
                         n 
                         ⁢ 
                         
                           p 
                           ^ 
                         
                         ⁢ 
                         
                           
                             F 
                             
                               
                                 2 
                                 ⁢ 
                                 n 
                                 ⁢ 
                                 
                                   p 
                                   ^ 
                                 
                               
                               , 
                               
                                 2 
                                 ⁡ 
                                 
                                   [ 
                                   
                                     
                                       n 
                                       ⁡ 
                                       
                                         ( 
                                         
                                           1 
                                           - 
                                           
                                             p 
                                             ^ 
                                           
                                         
                                         ) 
                                       
                                     
                                     + 
                                     1 
                                   
                                   ] 
                                 
                               
                             
                             
                               - 
                               1 
                             
                           
                           ⁡ 
                           
                             ( 
                             
                               a 
                               / 
                               2 
                             
                             ) 
                           
                         
                       
                     
                   
                   &lt; 
                   p 
                   &lt; 
                   
                     
                       
                         ( 
                         
                           
                             n 
                             ⁢ 
                             
                               p 
                               ^ 
                             
                           
                           + 
                           1 
                         
                         ) 
                       
                       ⁢ 
                       
                         
                           F 
                           
                             
                               2 
                               ⁡ 
                               
                                 [ 
                                 
                                   
                                     n 
                                     ⁢ 
                                     
                                       p 
                                       ^ 
                                     
                                   
                                   + 
                                   1 
                                 
                                 ] 
                               
                             
                             , 
                             
                               2 
                               ⁢ 
                               
                                 n 
                                 ⁡ 
                                 
                                   ( 
                                   
                                     1 
                                     - 
                                     
                                       p 
                                       ^ 
                                     
                                   
                                   ) 
                                 
                               
                             
                           
                           
                             - 
                             1 
                           
                         
                         ⁡ 
                         
                           ( 
                           
                             1 
                             - 
                             
                               a 
                               / 
                               2 
                             
                           
                           ) 
                         
                       
                     
                     
                       
                         n 
                         ⁡ 
                         
                           ( 
                           
                             1 
                             - 
                             
                               p 
                               ^ 
                             
                           
                           ) 
                         
                       
                       + 
                       
                         
                           ( 
                           
                             
                               n 
                               ⁢ 
                               
                                 p 
                                 ^ 
                               
                             
                             + 
                             1 
                           
                           ) 
                         
                         ⁢ 
                         
                           
                             F 
                             
                               
                                 2 
                                 ⁡ 
                                 
                                   [ 
                                   
                                     
                                       n 
                                       ⁢ 
                                       
                                         p 
                                         ^ 
                                       
                                     
                                     + 
                                     1 
                                   
                                   ] 
                                 
                               
                               , 
                               
                                 2 
                                 ⁢ 
                                 
                                   n 
                                   ⁡ 
                                   
                                     ( 
                                     
                                       1 
                                       - 
                                       
                                         p 
                                         ^ 
                                       
                                     
                                     ) 
                                   
                                 
                               
                             
                             
                               - 
                               1 
                             
                           
                           ⁡ 
                           
                             ( 
                             
                               1 
                               - 
                               
                                 a 
                                 / 
                                 2 
                               
                             
                             ) 
                           
                         
                       
                     
                   
                 
                 ) 
               
             
             = 
             
               1 
               - 
               
                 a 
                 . 
               
             
           
         
       
     
     Process  1200  further includes returning the critical value x and the confidence interval 1−α for p (step  1212 ). The critical value is found by taking the smallest value that will result in a fraction of samples less than x to be greater than p. x is used for static factor change detection by process  900  of  FIG. 9  and static factor change module  1100  of  FIG. 11 . For example, x is used as the critical value in step  902  of process  900 . The difference calculated through each iteration of steps  906 - 908  is compared to the value of x in step  910 . 
     Referring now to  FIG. 13 , a process  1300  of determining a P-value is shown, according to an exemplary embodiment. Process  1300  may be used to determine the P-value to use during baseline period detection of static factor changes (e.g., step  814  of  FIG. 8A ). In process  800  of  FIG. 8A , since the P-value is not known, an approximate P-value is found using process  1300 , and it is determined if the approximate P-value is suitable for use in process  800 . Compared to process  1200  for finding a critical value given a known P-value, process  1300  includes having a maximum of a statistic and determining the P-value instead. 
     Process  1300  includes determining an allowable tolerance on the approximation of the P-value (step  1302 ). The allowable tolerance may be set by a user of the building management system or may be determined by process  800  or another system of the building management system. The allowable tolerance relates to a confidence interval or region for which the P-value should fall under. For example, the set allowable tolerance may be compared to the width of the calculated confidence interval (see step  1308 ). 
     Process  1300  includes drawing a number of samples N of the distribution (e.g., data from the two data windows) at random (step  1304 ). The number of samples taken may be determined by approximation by taking the maximum of a normal random vector with appropriate correlations or by simulation of the actual process of calculating the maximum statistic. 
     Using the samples N, an approximate P-value is calculated as well as a confidence region or interval for the P-value (step  1306 ). The P-value is calculated using the method of step  814  of  FIG. 8A , according to an exemplary embodiment. The confidence region for the approximate P-value is also determined in step  1306 . The confidence region is a region for which the probability of the actual P-value being within the confidence region is greater than a specified threshold. The confidence region is a function of both the number of samples used in the Monte Carlo simulation (e.g., determined in step  1304 ) and the P-value (determined in step  1306 ). For example, for a given P-value P′ and confidence interval for the P-value of [P′−a,P′+b] (where a and b are chosen or calculated based on the desired threshold) the probability that the actual P-value falls within the confidence interval should be equal to a specified threshold (e.g., 95%). 
     Process  1300  further includes comparing the width of the confidence region to the tolerance determined in step  1302  (step  1308 ). If the confidence region is smaller than the tolerance, then the approximate P-value for which the confidence region is based on is accepted and the P-value is returned to, for example, process  800  for determining a change in static factors (step  1310 ). If the confidence region is greater than the tolerance, then the approximate P-value is rejected since the approximation of the P-value is not reliable or close enough. Process  1300  then draws another N samples for calculating a new approximate P-value (step  1304 ). Step  1308  of comparing the confidence region to the tolerance can be described as determining if enough samples N were used to approximate the P-value. For example, as more samples are used to approximate the P-value, the confidence region is reduced as the approximation is more accurate due to the use of additional samples. 
     One embodiment of the present disclosure relates to a method for calculating a P-value associated with a baseline period test statistic (e.g., a Z-statistic related to a difference in model coefficients). The method includes calculating an approximate P-value. The approximate P-value may be calculated using a number of random samples of the distribution of data used to calculate the model coefficients. The method further includes calculating a confidence interval for the P-value and determining the accuracy of the approximate P-value using the confidence interval. The method further includes estimating that the approximate P-value is a reliable P-value when the approximate P-value is within the confidence interval. The approximate P-value may then be used further steps and calculations (e.g., to determine energy usage and energy usage changes). If the approximate P-value is not within the confidence interval during the estimating step, then a new approximate P-value is calculated. The new approximate P-value may be calculated by selecting new samples of the distribution of data used to calculate the model coefficients. 
     One embodiment of the present disclosure includes a method for calculating a critical value associated with a known distribution of data points. The critical value is calculated such that the probability that any given test statistic is greater than the critical value matches a predefined probability of falsely determining that a model has changed in a significant way. The method includes determining a desired probability (e.g., user-selected, preselected). The method further includes simulating (e.g., using a Monte Carlo simulation), for varying samples of the known distribution of data relevant data points, a process of finding the test statistic. The results of the simulations, the desired probability, and a desired confidence interval are used to calculate the critical value. The method further includes using the critical value and confidence interval in a comparison against a test statistic to determine if the model has changed. 
     While many of the embodiments of the present disclosure relate to baseline models for energy use, measurement and verification features provided by the above systems and methods may be performed on many building resources or energy types. For example, the measurement and verification as described in the present disclosure may be applied to resource consumption types including electrical energy (or power/demand) usage, natural gas usage, steam usage, water usage, and other types of resource usage or consumption. As one example, gas and steam consumption can be used in the energy usage models described herein. As another example, water usage or consumption (and/or other types of resource consumption by a building) may be represented and processed using, e.g., the baseline models and varying related systems and methods of the present disclosure. 
     Using Uncertainty to Shorten the Measurement and Verification Process 
     In some circumstances, an ECM provider may specify the amount of savings (in energy or monetary terms) that is expected to result from retrofitting a building with energy conservation measures. Savings may be specified in per period terms (e.g., savings per year) or in other terms (e.g., savings over the life of the equipment, building, contract, repayment period). The terms may be specified in a contract between an ECM provider and an entity that is implementing ECMs within its building. For example, a contract may specify savings of 2 MWh per year relative to what would have been experienced with pre-retrofit equipment. In some cases, the ECM provider may be compensated based on an amount of savings relative to a previous system. 
     A period of monitoring and verification (M&amp;V) often follows a building retrofit. The M&amp;V process tracks the energy usage of a building after the retrofit (e.g., to check whether expected savings are realized). The duration of an M&amp;V period is usually agreed upon when the ECMs are implemented. For example, the length of the M&amp;V period (e.g., 6 years, 10 years) may be provided for in a contract. The M&amp;V period may be longer than necessary to demonstrate that the ECMs are providing the expected savings. 
     M&amp;V can be expensive, costing as much as 10% of the total project cost. The costs may be associated with a somewhat manual verification of the energy usage and savings over the M&amp;V period. In some cases, these costs are carried by the ECM provider rather than a customer. Ending a M&amp;V period sooner, when the savings are demonstrated to a particular degree of confidence, may advantageously minimize the costs (e.g., for the ECM provider, for the consumer of the ECM). 
     As discussed in above, the savings that may result from energy conservation measures can be calculated by subtracting the actual, measured energy consumption (with ECMs, after retrofitting) from the estimated baseline energy consumption (what would have occurred without ECMs). The estimated savings may be associated with a measure of uncertainty. That is, the estimated savings are known within a confidence interval. For example, the energy savings in the first month of the measurement and verification process may be 190,755 Wh±88,131 Wh with a confidence level of 95%. The uncertainty in the savings may arise from, e.g., the uncertainty in baseline model used to compute the savings. 
     Uncertainty and statistical confidence may be utilized to determine the end of an M&amp;V period. At the beginning of the M&amp;V period, the uncertainty (as a percentage of savings, baseline usage, etc.) in the calculated energy savings may be high. As the M&amp;V period continues, more and more data is obtained and the percent uncertainty should decrease. In some embodiments, once the contractually specified savings have been demonstrated to a particular degree of statistical confidence, the M&amp;V period may stop. In other embodiments, once the savings have been shown not to be realized, with a particular degree of statistical confidence, the M&amp;V period may stop. In some embodiments, the required confidence level for continuing M&amp;V and/or for discontinuing M&amp;V may be contractually specified. 
     Referring to  FIG. 14 , a block diagram of a building management system (BMS) computer system  200  for using uncertainty to shorten an M&amp;V process is shown, according to an exemplary embodiment. The BMS computer system  200  shown in  FIG. 14  may be similar to the BMS computer system described in the discussion of  FIG. 10  and like reference numbers are used for clarity. For example, BMS computer system  200  may contain a processing circuit  250 , including a processor  252  and memory  254 . Processor  252  and memory  254  may be similar to the processor and memory described in the discussion of  FIG. 10 . The BMS computer system  200  shown in  FIG. 14  includes additional components beyond those shown in  FIG. 10 . Those additional components will be described in the text to follow. 
     Memory  254  of  FIG. 14  is shown to include an uncertainty calculator  1400 . Uncertainty calculator  1400  is shown to receive information from baseline model  220  and data from energy savings data  226 . For example, uncertainty calculator  1400  may receive baseline energy usage data from baseline model  220  and energy savings per period from energy savings data  226 . In other embodiments, uncertainty calculator  1400  may additionally receive input from post retrofit variable data application  222  and measured energy consumption  224 . 
     Uncertainty calculator  1400  may be configured to calculate the uncertainty in the baseline model. Uncertainty is a mathematical quantity that describes the precision and accuracy of, e.g., the baseline model. In order to compute uncertainty, uncertainty calculator may account for, e.g., modeling errors, sampling errors, measurement errors, and errors relating to multiple distributions in the data. These potential sources of error are discussed in greater detail above. In some embodiments, uncertainty calculator  1400  may be configured to calculate the uncertainty in the actual, measured energy usage. Uncertainty calculator  1400  may also compute the uncertainty associated with energy savings. 
     Uncertainty calculator  1400  may be configured to compute statistical quantities, such as the confidence interval, associated with the baseline model and the energy savings. In some embodiments, uncertainty calculator  1400  may compute quantities related to uncertainty in the units of energy (e.g., Wh). In other embodiments, uncertainty calculator  1400  may be further configured to compute the uncertainty as a percentage of the baseline and/or as a percentage of the savings. In some embodiments, the statistical quantities relating to uncertainty may be calculated on a per period basis (e.g., per day, per month, etc.). Once calculated, uncertainty calculator  1400  may transmit the uncertainty quantities to uncertainty database  1402 . 
     In some embodiments, uncertainty calculator  1400  may be configured to compute a target savings per period during which ending M&amp;V is considered. Target savings may be considered by uncertainty calculator  1400  as spread out over a number of periods. For example, if the target energy savings is 2 MWh per year and ending M&amp;V is considered monthly, then the uncertainty calculator  1400  may utilize target savings of 
                 2   ⁢           ⁢   MWh       12   ⁢           ⁢   months       =     0.167   ⁢           ⁢   MWh   ⁢     /     ⁢   month           
in uncertainty calculations.
 
     Uncertainty calculator  1400  may also compute a cumulative target savings. Cumulative target savings is the sum of target savings to date (i.e., the present period and all previous periods). Continuing the example above, the cumulative target savings in the third month are
 
0.167MWh+0.167MWh+0.167MWh=0.5MWh
 
Uncertainty calculator  1400  may transmit target savings and cumulative target savings to uncertainty database  1402 .
 
     In some embodiments, savings, the associated uncertainty, target savings, and other quantities may be calculated at different intervals (e.g., daily, weekly, etc.). For example, savings, uncertainty, target savings, and other quantities may be calculated daily. Using the values in the example above, the daily target savings would be: 
                   2   ⁢           ⁢   MWh       365   ⁢           ⁢   days       =     0.005   ⁢           ⁢   MWh   ⁢     /     ⁢   day       ,         
and the cumulative target savings for the third month (e.g., March) would be:
 
(0.005MWh/day×31 days)+(0.005MWh/day×29 days)+(0.005MWh/day×31 days)=0.4986MWh.
 
As described below, these target savings values may variously be used in uncertainty calculations or in uncertainty-based determinations.
 
     Uncertainty calculator  1400  may be configured to compute upper confidence bounds (UCBs) and lower confidence bound (LCBs) associated with energy savings. UCBs and LCBs describe the margin of error associated with energy savings. That is, the energy savings for a given period is expected to fall within the UCB and LCB for the period, with a particular confidence level. UCBs may be defined as energy savings+uncertainty (with uncertainty expressed in terms of upper end energy savings margin for a given confidence level), and LCBs may be defined as energy savings−uncertainty (with uncertainty expressed in terms of lower end energy savings margin for a given confidence level). Uncertainty and energy savings may vary for each period they are calculated. Uncertainty calculator  1400  may transmit UCBs and LCBs to uncertainty database  1402 . UCBs and LCBs may be used in a variety of uncertainty-based calculations and determinations. 
     Memory  254  is shown to contain uncertainty database  1402 . Uncertainty database  1402  stores information associated with calculated uncertainty or variables that are used to calculate uncertainty. For example, uncertainty database  1402  may store a measure of uncertainty associated with the energy savings for a given period. Uncertainty database  1402  may be organized by monitoring period. For example, if uncertainty is calculated monthly, then uncertainty database  1402  may store variables such as monthly baseline energy usage, monthly actual energy usage, monthly energy savings, monthly uncertainty, monthly target savings, monthly UCB, monthly LCB, or other statistics used to convey or calculate uncertainty. 
     Memory  254  is shown to contain M&amp;V ending calculator  1404 . M&amp;V ending calculator  1404  is configured to determine the end point of an M&amp;V period based on the statistical quantities associated with uncertainty. Statistical quantities stored in uncertainty database  1402  are transmitted to M&amp;V ending calculator  1404 . M&amp;V ending calculator  1404  may also receive data from agreement information  232 . Such data could include the statistical confidence level for energy savings specified in the contract. 
     At the beginning of an M&amp;V period, a moving LCB may be less than cumulative target savings because the initial uncertainty in energy savings calculations is often relatively large (e.g., due to a small post ECM data set). In other words, initial energy savings is often associated with a large margin of error relative to the total accumulated savings. As the M&amp;V period progresses, the margin of error decreases as more and more data leads to smaller uncertainty. At some point (e.g., if savings are actually being realized, cumulative target savings was appropriately estimated at the beginning of a project, and uncertainty continues to decrease), the LCB should become greater than the cumulative target savings. When LCB is greater than cumulative target savings, the energy savings are known to be realized with a particular confidence level. The particular confidence level, e.g., 95%, may be contractually specified. In an example using 95%, the M&amp;V period may be ended because the probability that the actual savings is greater than the LCB is 0.95, and LCB is greater than the cumulative target savings. (LCB and cumulative target savings, in chart and graph form, are presented and discussed in the descriptions of  FIG. 18  and  FIG. 19B , below.) 
     In some embodiments (e.g., when the savings are not being realized), the UCB may become less than the cumulative target savings. An ECM provider may stop performing M&amp;V, prior to the end of an M&amp;V period, when it is statistically proven that the target savings are not being realized. This may advantageously minimize costs associated with continuing M&amp;V for an ECM provider. When the UCB is less than cumulative target savings, the energy savings are known not to be realized. When this is true, even the greatest amount of savings that could be realized is less than the target savings. 
     The M&amp;V ending calculator  1404  may also be configured to determine that the M&amp;V period should be continued, despite LCB being greater than cumulative target savings, if the M&amp;V period has not continued for a full cycle or for another minimum predetermined period of time. For example, a full cycle could be one full year or multiple full years. A full cycle may advantageously account for seasonal changes occurring during that time span, which may affect building energy usage differently. In other embodiments, a minimum period may be based on a specified number of energy days (cooling energy days or heating energy days), degree days (cooling degree days or heating degree days), or other parameters. Energy days and degree days are described in above sections with respect to previous Figures. The duration of a full cycle may be contractually specified, user specified, or dynamically adjusted by the system (e.g., in the event changes in the baseline model have been detected which may impact future performance). 
     The output of M&amp;V ending calculator  1404  (e.g., “continue M&amp;V,” “end M&amp;V,” etc.) may be transmitted to facility monitoring module  230 . Facility monitoring module  230  may be similar to the facility monitoring module described in the discussions of  FIGS. 10 and 11 , above. For example, facility monitoring module  230  may transmit and receive data to/from user output/system feedback module  234 . User output/system feedback module  234  may be configured to output an audio and/or visual indication or alert (e.g., on an electronic display) to at least one of a memory device, a user device, or another device on the building management system. The audio and/or visual indication may reflect the status of the M&amp;V process as determined by M&amp;V ending calculator  1404 . M&amp;V ending calculator may also transmit its output to uncertainty database  1402 . 
     In other embodiments, user interfaces (e.g., graphical user interfaces) generated by user output system feedback module  234  are configured to display uncertainty calculations, related graphs or tables (e.g., as shown in subsequent Figures), or provide other dashboards. One such dashboard may use the output of M&amp;V calculator  1404  to provide a suggestion to a user (e.g., a building manager, a contract manager) to end M&amp;V. The dashboard may include options for allowing the user to accept or deny the suggestion. Moreover, the dashboard may allow the user to specify an extension to the M&amp;V system (e.g. M&amp;V to run one more year under a condition where the M&amp;V ending calculator  1404  suggests ending). 
     Referring to  FIG. 15 , a baseline calculation module  210  for calculating uncertainty is shown in greater detail, according to an exemplary embodiment. Baseline calculation module  210  may be a similar to the baseline calculation module described in the discussion of  FIG. 2 . Baseline calculation module  210  of  FIG. 15  may contain a data clean-up module  212 , data synchronization module  214 , regression analysis module  216 , and cross-validation module  218 . These components may be similar to the corresponding components described in the discussion of  FIG. 2 . 
     Baseline calculation module  210  may include uncertainty module  1500  (also shown in  FIG. 14 ). Uncertainty module  1500  may be configured to collect and store data relating to potential sources of error in calculating the baseline model. Types of error may include modeling errors, sampling errors, measurement errors, and errors relating to multiple distributions in the data. Uncertainty module  1500  may transmit data relating to error first to baseline model  220  and then to uncertainty calculator  1400  ( FIG. 14 ). The data may be transmitted in a format that allows for uncertainty calculator  1400  to estimate the uncertainty relating to the baseline model generation. 
     Referring to  FIG. 16 , a flow diagram of a process  1600  for using uncertainty to shorten the M&amp;V process is shown, according to an exemplary embodiment. Process  1600  may be implemented on, e.g., processing circuit  250  of BMS computer system  200  ( FIG. 14 ). The calculations required by process  1600  may be completed by uncertainty calculator  1400  and/or M&amp;V ending calculator  1404  and may be stored in uncertainty database  1402 . User output/system feedback module  234  may transmit the end result of process  1600  to a user. The results of process  1600  may also or alternatively be communicated to another system via communications interface  242 . 
     Process  1600  includes setting a confidence level ( 1602 ). Confidence level is a statistical description of the probability that an estimated value exists within a certain range. A quantity, such as energy savings resulting from ECMs in a building, may be known to be within a given range with specific confidence level. For example, energy savings in a particular month may be 190,755 Wh±88,131 Wh with a confidence level of 95%. The confidence level may indicate that, if the data used to calculate energy savings for the particular month were collected again, the energy savings would fall within the specified margin of error 95% of the time. Various confidence levels may be chosen (e.g., 68%, 75%, 95%, 99%, etc.). A higher confidence level may require a greater amount of data (i.e., a longer M&amp;V period) to verify that energy savings are being realized (and vice versa). The confidence level may be specified in a contract, selected based on a budget, adjusted over time, user selected, or otherwise determined. In various embodiments, the confidence levels for ending M&amp;V after proving successful savings and ending M&amp;V after proving a savings failure are the same or different. 
     Process  1600  includes steps  1604 - 1610  that are repeated for each period in which energy savings and the related uncertainty are calculated. A user may specify how frequently energy savings are calculated (e.g., every day, every month, etc.). 
     Process  1600  includes calculating cumulative target savings ( 1604 ). Cumulative target savings may be described as a target amount of energy savings resulting from energy conservation measures in a building. In some embodiments, a user interface for allowing an operator of the building management system to set or adjust the target amount may be provided. 
     Process  1600  includes calculating the lower confidence bound (LCB) (step  1606 ). Lower confidence bound may be described as a least amount of energy savings (for a given confidence level) resulting from energy conservation measures in a building. The LCB for a given period is defined as energy savings for that period−uncertainty for that period. In other words, the lower confidence bound describes, with a given confidence value, the lowest value the energy savings may be were the same data collected again. The LCB will likely be different for each period because uncertainty and energy savings are usually both different for each period. In some embodiments, the upper confidence bound (UCB) is calculated. Upper confidence bound may be described as a greatest amount of energy savings (for a given confidence level) resulting from energy conservation measures. The UCB for a given period is defined as energy savings for that period+uncertainty for that period. 
     Process  1600  includes comparing the LCB to the cumulative target savings (step  1608 ). Comparing the LCB to the cumulative target savings may indicate whether or not the contractually specified (or otherwise target) savings are being realized at a particular confidence level. The comparison may be a straight comparison completed using the LCB and cumulative target savings for a particular period. In other embodiments, the comparison may include or be preceded by a transformation of one or both values. In some embodiments, the UCB is compared to the cumulative target savings. This comparison may indicate whether the target savings are not being realized at a particular confidence level. 
     Process  1600  includes continuing the M&amp;V period when the LCB is less than cumulative target savings (step  1610 ). If the LCB is less than cumulative target savings, the lowest amount the energy savings could be for that period (with the specified statistical confidence) is less than the cumulative target savings. Thus, the actual energy savings may be less than the contractually specified savings. The M&amp;V process can be continued to include additional periods of data to determine whether the cumulative target savings condition is satisfied in the future. In some embodiments, the M&amp;V period is continued when the UCB is greater than the cumulative target savings. 
     Process  1600  include ending the M&amp;V period when the LCB is greater than cumulative target savings (step  1612 ). If the LCB is greater than the cumulative target savings, then even the lowest the energy savings could be for that period (with the specified statistical confidence) is greater than the cumulative target savings. When the LCB is greater than the cumulative target savings, the actual energy savings to date may be much greater than the cumulative target savings. Thus, the contractual energy savings condition likely has been satisfied. In some embodiments, the M&amp;V period is ended when the UCB is less than cumulative target savings. If the UCB is less than the cumulative target savings, then even the greatest the energy savings could be for that period (with the specified statistical confidence) is less than the cumulative target savings. Thus, the contractual energy savings condition likely has not been satisfied. 
     Advantageously, the absolute amount of uncertainty may not be critical when comparing LCB and cumulative target savings. More important may be the fact that the uncertainty in energy savings is known with a given confidence level (set in step  1602 ). Thus, when the LCB is greater than cumulative target savings, the energy savings are known to satisfy the contractual energy savings condition with a particular confidence level. Similarly, when the UCB is less than the cumulative target savings, the energy savings are known not to satisfy the contractual energy savings condition with a particular confidence level. 
     Process  1600  includes alerting the user that the M&amp;V period may end (e.g., because the target savings are being realized, because the target savings are not being realized) ( 1614 ). The BMS may be configured to output an indication that the M&amp;V period has ended (e.g., or is recommended to end, may end, etc.) to at least one of a memory device, a user device, or another device on the building management system. The indication may be an audio or visual representation on a BMS user interface. In some embodiments, ending the measurement and verification period does not require user input of the ending to be received from an operator of the building management system. In other embodiments, the BMS may present the user with the choice of ending M&amp;V, continuing M&amp;V, or otherwise adjusting the M&amp;V process. 
     Referring to  FIG. 17 , a more detailed flow diagram of a process  1700  for using uncertainty to shorten the M&amp;V process is shown, according to an exemplary embodiment. Process  1700 , like process  1600 , may be implemented on, e.g., processing circuit  250  of BMS computer system  200  ( FIG. 14 ). Uncertainty calculator  1400  and/or M&amp;V ending calculator  1404  conduct calculations and store the results in uncertainty database  1402 . The end result of process  1700  may be output to a user via user output/system feedback module  234 . 
     Process  1700  includes setting the target savings ( 1702 ). The target savings may be those predicted to result from a retrofit of a building with EOMs. The target savings may be specified in a contract. According to some embodiments, the savings may be described in per period terms (e.g., savings per year). For example, a contract may specify savings of 2 MWh per year. Process  1700  includes setting the monitoring period length ( 1704 ). The length of the monitoring period describes how frequently the energy savings, uncertainty, and related quantities (uncertainty as a percentage of baseline, uncertainty as a percentage of savings, LCB, UCB, etc.) are calculated. For example, energy savings could be calculated daily. In some embodiments, the intervals of time when energy savings are calculated may be different from the intervals when ending the M&amp;V process is considered. For example, energy savings may be calculated daily, while ending the M&amp;V period may be done monthly. In other embodiments, energy savings and ending the M&amp;V process may be considered at the same intervals. Process  1700  further includes calculating the target savings per monitoring period ( 1706 ). The target savings per monitoring period may be calculated by equally dividing the target savings (e.g., 2 MWh/year) over the monitoring periods (e.g., 365 days). 
     Referring still to  FIG. 17 , process  1700  is shown to include setting confidence level ( 1708 ). The confidence level is a statistical indicator of certainty, as described in the discussion of  FIG. 16 , above. 
     In process  1700 , steps  1710 - 1726  are repeated for each period. A period may describe the intervals between which energy savings and related quantities are calculated, and/or the intervals between which ending the M&amp;V process is considered. 
     Process  1700  includes calculating the cumulative target savings ( 1710 ). Cumulative target savings for a particular period is the sum of the target savings (calculated in step  1706 ) for that period and all past periods (with period length set in step  1704 ). Cumulative target savings may be compared to the energy savings lower confidence bound (in step  1722 ) as a part of a determination of whether the M&amp;V process can be ended. 
     Process  1700  is shown to include steps  1712 - 1716  to calculate the energy savings. Steps  1712 - 1716  may be similar to those described in the discussion of  FIGS. 1A-1B , above. 
     Process  1700  further includes calculating a mathematical uncertainty in the energy savings ( 1718 ). Uncertainty in the energy savings may be the result of errors in the baseline model used to calculate energy savings. Uncertainty may be calculated as described in the discussion  FIG. 14 , above. 
     Process  1700  includes calculating a lower confidence bound ( 1720 ). The LCB may be calculated as described in the discussion of  FIG. 16 , above. As noted above, process  1700  includes comparing lower confidence bound to cumulative target savings ( 1722 ). This comparison may indicate whether target energy savings is being realized to a particular level of confidence. In some embodiments, the upper confidence bound may be calculated and compared to the cumulative target savings. 
     Process  1700  includes continuing the M&amp;V process when the lower confidence bound is less than cumulative target savings ( 1726 ). As described in the discussion of  FIG. 16 , above, actual energy savings may be less than the expected savings (set in step  1702 ) when LCB is less than the cumulative target savings. The M&amp;V period may be continued to determine if the energy savings will meet or exceed the target savings in the future. Similarly, the M&amp;V process may be continued when the upper confidence bound is greater than the cumulative target savings. 
     Process  1700  includes determining whether the M&amp;V process has completed a full cycle when the lower confidence bound is greater than cumulative target savings ( 1724 ). As described in the discussion of  FIG. 16 , above, if the LCB is greater than the cumulative target savings, then the contractual energy savings condition may be satisfied, to a particular level of statistical confidence (set in step  1708 ). The LCB may become greater than cumulative target savings before M&amp;V process has completed a full cycle. A full cycle of a monitoring period may be based on, e.g., time, number of energy days, degree days, etc., as described in the discussion of  FIG. 14 , above. 
     Process  1700  includes continuing the M&amp;V process ( 1726 ) if the M&amp;V process has not completed a full cycle and ending the M&amp;V process ( 1728 ) if the M&amp;V process has completed a full cycle. For example, for an M&amp;V period that begins in January, the lower confidence bound may be greater than cumulative target savings by August of the same year. Because a full cycle (in this case, a full year) has not been completed, the M&amp;V period can be continued through January of the following year. If, for example, lower confidence bound becomes greater than cumulative target savings after the M&amp;V process has completed a full cycle, the M&amp;V process can end. 
     Process  1700  includes alerting the user when the target savings are being realized and a full M&amp;V cycle has been completed ( 1730 ). The user may be alerted as described in the discussion of  FIG. 16 , above. For example, a visual indication may be presented on a user device. According to some embodiments, user input may be required to continue or end the M&amp;V process. In yet other embodiments, the M&amp;V process may end automatically. 
     Referring to  FIG. 18 , a chart of energy savings and statistical quantities related to uncertainty over a period of time is shown. The chart of  FIG. 18  represents a set of data collected over the M&amp;V period. In some embodiments, the end of an M&amp;V period (and/or other M&amp;V decisions) may be determined using a database set or chart such as that of  FIG. 18 . The chart (i.e., database, array, data object, etc.) may be created, maintained, updated, processed, and/or stored in uncertainty database  1402  ( FIG. 14 ). The calculations reflected on the chart of  FIG. 18  may be completed by uncertainty calculator  1400  and/or M&amp;V ending calculator  1404  ( FIG. 14 ). 
     The chart of  FIG. 18  may be output to at least one of a memory device, a user device, or another device on the building management system. The output may be a graphical user interface (e.g., on a client device, on a mobile device, generated by a web server, etc.). For example, the chart of  FIG. 18  may be output via user output/system feedback module  234  ( FIG. 14 ). Depending on the embodiment, the chart may include a portion or all of the columns and rows of data shown in  FIG. 18 . In some embodiments, a user may be able, via a user interface, to choose which columns and rows should be visible. A user may also be able to add columns and/or rows of data not shown in  FIG. 18 , including, e.g., weather, operating condition of the building (e.g., occupied or unoccupied), or other parameters. In other embodiments, the process will operate without displaying a graphical representation of a chart. While the chart of  FIG. 18  is shown as a two dimensional chart, another information structure suitable for representing and storing the data of the chart may be used. For example, a relational database having one or more related tables may be used. 
     The chart of  FIG. 18  include fields populated by dates for period start  1800  and period end  1802 . The length of a period (e.g., day, week, month) for which energy savings and the associated quantities are calculated may vary depending on user selection, contract terms, data collection frequency or other dependencies. In the illustration of  FIG. 18 , each row represents an additional month of data and the rows contain cumulative data. That is, each period start begins on the same date (January 1) and ends on different dates (e.g., February 1, March 1, etc.). In other embodiments, data may be presented on a “per period” (and not cumulative) basis. 
     The chart of  FIG. 18  is shown to include baseline usage  1804 . Baseline usage  1804  may be the estimated energy usage of the building without ECMs. In the embodiment of  FIG. 18 , baseline usage is displayed in watt hours (Wh). The baseline usage is predicted from the baseline model, as described above. The chart of  FIG. 18  also includes actual energy usage  1806 . Actual energy usage is the measured energy usage of the building with ECMs. In the embodiment of  FIG. 18 , actual energy usage is displayed in watt hours (Wh). Other units may be used in different embodiments. 
     The chart of  FIG. 18  also includes energy savings  1808 . Energy savings are the difference between energy usage predicted by the baseline model  1804  and the actual, measured energy usage  1806 . In the embodiment of  FIG. 18 , energy savings are displayed in watt hours (Wh). 
     The chart of  FIG. 18  includes savings as a percentage of baseline energy usage  1810 . Savings as a percentage of baseline  1810  may represent how much energy was actually expended by the building with the ECMs versus how much energy would have been consumed by the building without ECMs. It should be appreciated that 0.15 represents, e.g., 15%. 
     The chart of  FIG. 18  includes uncertainty  1812 . As described above, uncertainty  1812  may be calculated by uncertainty calculator  1400  ( FIG. 14 ) and calculated with reference to a given confidence level. In the embodiment of  FIG. 18 , the confidence level used was 95%. 
     Uncertainty  1812  may be the result of various errors in baseline usage  1804  or actual usage  1806 . In some embodiments, the uncertainty as result of errors in the baseline model and uncertainty as a result of errors in the actual usage may be displayed separately (e.g., as different columns). In the embodiment of  FIG. 18 , uncertainty is displayed in watt hours (Wh). Other units may be used in different embodiments. 
     The chart of  FIG. 18  includes uncertainty as a percentage of the baseline usage  1814  and as a percentage of savings  1816 . These percentages may reflect a decrease in uncertainty percentages as more data is included in calculating uncertainty (i.e., when the M&amp;V period continues). 
     The chart of  FIG. 18  includes the number of cooling energy days 1818. Cooling energy days are, as described above, a predictor variable that combines both outside temperature and outside air humidity. As reflected in the chart, the number of cooling energy days, per month, may be greatest during the summer months. The chart may also or alternatively include, e.g., the number of heating energy days, cooling degree days, heating degree days, etc. 
     The chart of  FIG. 18  includes the cumulative target savings  1820 . As discussed above, the cumulative target savings are the sum of target savings up to and including the month in which the energy savings are calculated. Target savings may be calculated by dividing a yearly target savings over the intervals when energy savings are calculated (e.g., monthly, weekly, daily). In the embodiment of  FIG. 18 , the chart reflects a yearly target savings of 2 MWh, with energy savings calculated daily. 
     The chart of  FIG. 18  includes lower confidence bound  1822  and upper confidence bound  1824 . Lower confidence bound  1822  may be the difference of savings  1808  and uncertainty  1812 . Upper confidence bound may be the sum of savings  1808  and uncertainty  1812 . Lower confidence bound  1822  and upper confidence bound  1824  may define the expected savings  1808 , with a given confidence level (95%, in the embodiment of  FIG. 18 ). 
     The chart of  FIG. 18  includes row  1826  illustrating an identified data set where it may be determined that the M&amp;V period can be ended. In row  1826 , the value of the lower confidence bound  1822  is greater than the cumulative target savings  1820  for the first time during the M&amp;V process. The data reflects that energy savings have been demonstrated beyond a particular confidence level (95%, in the embodiment of  FIG. 18 ). In some embodiments, the M&amp;V period may end after August 1 (the time period of the data in row  1826 ). In other embodiments, the M&amp;V period may continue until a full cycle (e.g., a full year) of data has been collected and the lower confidence bound has been maintained above the cumulative target savings. When energy savings are not being realized, the chart of  FIG. 18  may shown that the value of the upper confidence bound becomes less than the cumulative target savings. 
     Referring to  FIGS. 19A-19C , plots of energy savings and related quantities, over time, are shown, according to exemplary embodiments. As described in the discussion of  FIG. 18 , the plots of  FIGS. 19A-190  may be created by uncertainty calculator  1400  and/or M&amp;V ending calculator  1404 , and stored in uncertainty database  1402  ( FIG. 14 ). The plots may be output on a graphical user interface (e.g., on a client device). In some embodiments, the plots may be output with the chart of  FIG. 18 . In various embodiments, the data in  FIGS. 19A-19C  may be combined into a single plot or may be spread across multiple plots. A user may specify the manner in which the data is organized and displayed on a graphical user interface. 
     Referring to  FIG. 19A , a plot of energy savings, LCB, UCB, and uncertainty as a percentage of energy savings, over time, is shown, according to an exemplary embodiment. Note that the plot of  FIG. 19A  contains two y-axes. The y-axis on the left (energy) corresponds to the scale for savings  1900 , LCB  1902 , and UCB  1904 . The y-axis on the right (uncertainty) corresponds to the scale for uncertainty as a percentage (e.g., where 0.15 uncertainty equals 15%) of savings  1906 . 
     The plot of  FIG. 19A  includes energy savings  1900 . Energy savings  1900  are the difference between baseline energy usage and actual energy usage.  FIG. 19A  also includes LCB  1902  and UCB  1904 . LCB  1902  is defined as energy savings−uncertainty, and UCB  1904  is defined as energy savings+uncertainty. In the example of  FIG. 19A , energy savings  1900 , LCB  1902 , and UCB  1904  are cumulative (from the beginning of the M&amp;V process). In other embodiments, energy savings  1900  may be presented on a per period (e.g., per month) basis. 
     The plot of  FIG. 19A  includes uncertainty as a percentage of savings  1906 .  FIG. 19A  reflects that, in the underlying data, uncertainty  1906  decreases as the M&amp;V process continues and more data is available to calculate uncertainty. While in  FIG. 19A  uncertainty is shown as a percentage of savings, in other embodiments the plot may illustrate uncertainty as a percentage of baseline usage. 
     Referring to  FIG. 19B , a plot of cumulative target savings  1908  and lower confidence bound  1902 , over time, is shown, according to an exemplary embodiment. Cumulative target savings  1908  are, as discussed above, the sum of target savings per monitoring period to date. Cumulative target savings  1908  and lower confidence bound  1902  may be compared to determine whether target savings are being realized (as described in  FIGS. 16 and 17 ). As reflected in  FIG. 19B , lower confidence bound  1902  is less than the cumulative target savings  1908  at the beginning of the M&amp;V process. This is because the uncertainty at the beginning of the M&amp;V process may be high, causing a wide margin (defined by the upper and lower confidence bounds) for energy savings. As the uncertainty decreases, the margin for energy savings becomes narrower (the upper confidence bound decreases and the lower confidence bound increases). When the lower confidence bound  1902  becomes equal to and greater than target savings (as shown near the middle of the plot), the target savings are known to be realized to a degree of statistical confidence. In the embodiment of  FIG. 19B , the confidence level may be 95%. In some embodiments, the plot of  FIG. 19B  may show the cumulative target savings and the upper confidence bound, over time. At the beginning of the M&amp;V process, the upper confidence bound may be greater than the cumulative target savings. When the energy savings are not being realized, the upper confidence bound becomes equal to and then less than the cumulative target savings. 
     Referring to  FIG. 19C , a plot of actual energy savings and cumulative target energy savings, over time, is shown, according to an exemplary embodiment. Actual energy savings  1900  is defined as baseline energy usage−measured energy usage. The plot reflects that actual savings  1900  were greater than cumulative target savings  1908  from the beginning of the M&amp;V process. However, the savings were only proven to a degree of statistical confidence when the lower confidence bound of the actual savings become greater than the cumulative target savings (as shown in  FIG. 19B ). When the target savings are not being realized,  FIG. 19C  may show that the actual energy savings are less than the cumulative target savings for at least a portion of the M&amp;V process. 
     In an exemplary embodiment, the plots of  FIGS. 19A-190  may be manually evaluated on a graphical user interface to make decisions with respect to uncertainty of M&amp;V and/or ending or otherwise revising an M&amp;V relationship. In other embodiments, the plots of  FIGS. 19A-190  may be evaluated by a computer module of the BMS. For example, an automated module may evaluate certain relationships between plots on the graphs. In one example, the crossing point of plots  1904 ,  1902 , or  1900  with respect to plot  1906  may be used to make an automated decision or recommendation. The system, for example, may determine and then state “uncertainty as a percent of savings has crossed savings and its uncertainty bands.” Such a determination may be used to end or shorten an M&amp;V process. In varying embodiments, such a determination (as well as the other uncertainty determinations described herein) may be used to monitor, gauge, or adjust contract performance in other ways. For example, when savings are sufficiently high and uncertainty is sufficiently low, one or more incentive thresholds may be hit. On the other hand, using the uncertainty-based calculations and determinations described herein a company offering EOMs may be able to provide potential customers with guarantees, claw-backs (in the event guarantees are not hit), or other incentives. 
     Configurations of Various Exemplary Embodiments 
     The construction and arrangement of the systems and methods as shown in the various exemplary embodiments are illustrative only. Although only a few embodiments have been described in detail in this disclosure, many modifications are possible (e.g., variations in sizes, dimensions, structures, shapes and proportions of the various elements, values of parameters, mounting arrangements, use of materials, colors, orientations, etc.). For example, the position of elements may be reversed or otherwise varied and the nature or number of discrete elements or positions may be altered or varied. Accordingly, all such modifications are intended to be included within the scope of the present disclosure. The order or sequence of any process or method steps may be varied or re-sequenced according to alternative embodiments. Other substitutions, modifications, changes, and omissions may be made in the design, operating conditions and arrangement of the exemplary embodiments without departing from the scope of the present disclosure. 
     The present disclosure contemplates methods, systems and program products on any machine-readable media for accomplishing various operations. The embodiments of the present disclosure may be implemented using existing computer processors, or by a special purpose computer processor for an appropriate system, incorporated for this or another purpose, or by a hardwired system. Embodiments within the scope of the present disclosure include program products comprising machine-readable media for carrying or having machine-executable instructions or data structures stored thereon. Such machine-readable media can be any available media that can be accessed by a general purpose or special purpose computer or other machine with a processor. By way of example, such machine-readable media can include RAM, ROM, EPROM, EEPROM, CD-ROM or other optical disk storage, magnetic disk storage or other magnetic storage devices, or any other medium which can be used to carry or store desired program code in the form of machine-executable instructions or data structures and which can be accessed by a general purpose or special purpose computer or other machine with a processor. Combinations of the above are also included within the scope of machine-readable media. Machine-executable instructions include, for example, instructions and data which cause a general purpose computer, special purpose computer, or special purpose processing machines to perform a certain function or group of functions. 
     Although the figures may show a specific order of method steps, the order of the steps may differ from what is depicted. Also two or more steps may be performed concurrently or with partial concurrence. Such variation will depend on the software and hardware systems chosen and on designer choice. All such variations are within the scope of the disclosure. Likewise, software implementations could be accomplished with standard programming techniques with rule based logic and other logic to accomplish the various connection steps, processing steps, comparison steps and decision steps.