Nonintrusive Load Monitoring of Variable Speed Drive Cooling Systems

To improve the energy efficiencies of building cooling systems, manufacturers are increasingly utilizing variable speed drive (VSD) motors in system components, e.g. compressors and condensers. While these technologies can provide significant energy savings, these benefits are only realized if these components operate as intended and under proper control. Undetected faults can foil efficiency gains. As such, it’s imperative to monitor cooling system performance to both identify faulty conditions and to properly inform building or multi-building models used for predictive control and energy management. This paper presents nonintrusive load monitoring (NILM) based “mapping” techniques for tracking the performance of a building’s central air conditioning from smart electrical meter or energy monitor data. Using a multivariate linear model, a first mapping disaggregates the air conditioner’s power draw from that of the total building by exploiting the correlations between the building’s line-current harmonics and the power consumption of the air conditioner’s VSD motors. A second mapping then estimates the air conditioner’s heat rejection performance using as inputs the estimated power draw of the first mapping, the building’s zonal temperature, and the outside environmental temperature. The usefulness of these mapping techniques are demonstrated using data collected from a research facility building on the Masdar City Campus of Khalifa University. The mapping techniques combine to provide accurate estimates of the building’s air conditioning performance when operating under normal conditions. These estimates could thus be used as feedback in building energy management controllers and can provide a performance baseline for detection of air conditioner underperformance.


I. INTRODUCTION
In developed countries, energy expenditures in buildings represent 20-40% of the total primary energy use [1]- [3]. As this energy typically derives from fossil fuels, buildings represent major emitters of carbon dioxide and leading contributors to climate change [3], [4]. On average, 19% of a building's energy expenditure goes to cooling and ventilation [5]. In hot The associate editor coordinating the review of this manuscript and approving it for publication was Ming Luo . and humid climates, e.g., the United Arab Emirates (UAE), cooling loads can exceed 60% on peak summer days [6].
A variety of potential solutions have been explored to improve cooling efficiencies in buildings including integrating variable speed/frequency drive (VSD/VFD) motors into system components [7]- [9]. These VSDs allow the cooling system, e.g., a centralized air conditioner or circulating chilled water system, to adjust its operational capacity to match the building's cooling needs [7], [10]. Further, with sufficient thermal storage and the cooling system's input for urban climate modeling. These results also show the potential for these mapping techniques to provide ''continuous commissioning'' of the cooling system. For example, the estimated heat rejection could be used as an input for building heat transfer models to generate an expected indoor temperature. If the true indoor temperature is not within an acceptable range around this estimate, that could indicate a cooling system fault, e.g. a refrigerant leak. This particular fault was observed during testing, and the estimated heat rejection during the fault is compared against the measured heat rejection to illustrate this point.

II. NONINTRUSIVE LOAD MONITORING REVIEW
The concept of nonintrusively detecting the operation of individual loads in a building by analyzing the characteristics of its utility power flow was first introduced by MIT researchers George Hart, Ed Kern, and Fred Schweppe in the 1980s [31]. In his seminal paper introducing nonintrusive load monitoring to the greater engineering community [32], Hart effectively provided a research roadmap for NILM development including outlining several overarching strategies for load disaggregation depending on the types of electric loads present and distinguishing between ''manualsetup'' and ''automatic-setup'' NILM systems. These setup categories forecast the use of supervised and unsupervised machine learning techniques, respectively, as Hart describes manual-setups as those requiring observing and labeling recorded data while turning loads on and off, and automatic-setups as those using a priori information for identifying loads.
Since [32], numerous teams have advanced the state of the art in NILM research, and several review papers exist to provide an overview of these advances, e.g. [22], [23], [33]- [36]. In particular, Zoha et al. [22] provided a very comprehensive review of NILM research at the time of its authorship, covering the various NILM frameworks used by researchers, the features extracted from current, voltage, and/or power measurements for load identification, and the supervised and unsupervised algorithms used to identify load operations. More recently, Bonfligli et al. [35] provided a review focused exclusively on unsupervised learning based NILM techniques and Pereira and Nunes [36] provided a review of publicly available datasets and toolkits for testing NILM techniques and the various metrics used to evaluating load identification efficacy.
Building loads generally fall into one of four categories: on/off loads, multi-state loads, constant base loads (always on), and continuously variable loads [22], [32], [37]. Historically, on/off loads (e.g. non-dimmable lights), multi-state loads (e.g. washing machines), and constant base loads (e.g. unswitched emergency lighting) have been the dominant categories. The first two types of loads produce characteristic transients and steady-state step change events in electrical measurements which can be used as ''signatures'' for detecting and identifying load operations [24]. Event-based techniques that utilize these signatures generally feature supervised learning load identification algorithms. These algorithms typically outperform unsupervised approaches in terms of accuracy, but require labelled data for training, which can be difficult or expensive to obtain. Numerous supervised learning algorithms have been applied for NILM, including support vector machines (SVMs) [38], [39], and both shallow [40], [41] and deep neural networks (NNs) [42], [43]. Researchers have also used pre-filtering techniques to generate input features for NNs across vastly different time-scales [37] and combined NN-based load identification with unsupervised learning optimization techniques to improve results [44].
Non-event based, unsupervised learning methods have the advantage of not requiring extensive labelled data and thus are more broadly applicable. Further, unsupervised approaches that analyze steady-state electrical measurements and use optimization techniques to infer the best-fit combination of loads to match the measurements can also account for constant base loads, e.g. [45], which uses an evolutionary optimization algorithm. Suzuki et al. used a different optimization technique incorporating integer programming to determine the active loads from the bulk current measurements [46]. Bhotto et al. later improved upon this work by incorporating state diagrams to correct integer programming algorithm outputs [47]. Another widely used unsupervised approach incorporates Factorial Hidden Markov Models (FHMMs), e.g. [48]- [50], to model the state transitions of various loads in a household or facility.
The final load category, continuously variable loads, are increasing in prevalence particularly as VSD motor loads become more cost effective. These types of loads do not draw power at only discrete levels, and thus they are generally not identifiable by either event-based techniques or non-event based techniques developed for identifying the other three load types. However, it is well understood that these sort of variable loads impart harmonic content into building line currents due to their front-end rectifiers. Previously, Lee et al. reported a correlation between the power draw of VSD driven blower fans in a commercial building and the 5th harmonic magnitude imparted on a single-phase current [51]. The authors used a piecewise power function model to estimate the fans' power draw from this 5th harmonic measurement. Wichakool et al. developed an analytical switching-function-based approach to estimate the power draw of VSDs and other rectified loads from three-phase current harmonics [52]. Wichakool later followed this with an empirical approach exploiting structural features in rectified current waveforms [53], however both of these techniques were only tested in lab-based setups with limited total loads.
In this paper, we build upon [51]- [53] by taking an empirical, machine learning-based approach to disaggregating the power consumption of multiple VSD cooling system loads in a commercial building. In this approach, the in-phase, quadrature, and apparent, odd-order current harmonics up to the 7th harmonic on all three building phases are measured and used as feature inputs to a multivariate linear estimator of total power draw for the VSD loads. Fig. 1 illustrates the basic operation of a typical vapor-compression cooling system like the one at the experimental site. In this system, refrigerant in a mixed liquid/vapor state passing through coils in an evaporator absorbs heat from within a facility (q c ) as warm indoor air (or water in the case of a chilled water heat exchanger) passes over the coils. This exchange of heat causes the refrigerant to change phases from a liquid to a gas (saturated vapor). The gaseous refrigerant then flows to a compressor unit outside the facility where it's compressed and coincidentally increases in temperature. The hot high-pressure vapor then enters a condenser which removes heat from the refrigerant vapor by blowing ambient air (or circulating water in the case of a water-cooled condenser) across the coils or tubes containing the refrigerant vapor. This causes the refrigerant to condense into a saturated liquid as heat (q h ) transfers from the refrigerant to the air or water. Finally, the liquid refrigerant passes through an expansion valve where it experiences an abrupt decrease in pressure and evaporates to a cool refrigerant liquid/vapor mixture, thus completing the cycle. The result is facility cooling by the evaporator with the absorbed heat transported via the refrigerant to the condenser where it's rejected to the environment along with the total work (W ) done on the system, i.e.,

III. BUILDING COOLING SYSTEMS
Taking the time-derivative of this equation restates the relationship in terms of the rates of heat transfer and the electric power delivered to the cooling system components (P in = dW dt ), i.e.,q

VARIABLE SPEED DRIVE COOLING SYSTEMS
The compressor and both the evaporator and condenser fans (or pumps in the case of water systems) are driven by electric motors. In conventional air conditioning systems, these motors are typically induction machines connected directly to the electrical supply. While simple and robust, these systems suffer from inefficiencies particularly during times of light cooling [10]. Inverter-driven, variable speed drive (VSD) motors can significantly improve cooling efficiencies [10], [12], [54], [55], particularly if integrated into the compressor as it draws the bulk of the system's electrical power. VSDs allow the motors to vary their speed, and therefore the system's cooling capacity. This enables low-lift operation, where the ratio of refrigerant condensing pressure to evaporating pressure is kept low, greatly reducing the required work while still providing a similar cooling effect [54], [55]. Inverter-driven VSD components typically employ DC-rectification of the utility line voltages as depicted in Fig. 2(a). An inverter then converts this DC voltage back to AC with the frequency of the voltage waveforms controlled via modulation of the inverter's power electronic switches [56]. As the motor draws current from the DC-link, nonlinear currents containing harmonic components are drawn from the AC lines [57] if active power factor correction is not deployed. This is the case at our experiment site. However, even systems that employ active wave-shaping create smaller but typically detectable harmonic signatures that can be used for tracking. Single-phase rectifiers impart a single positive peak and single negative peak onto the line currents during each line cycle ( Fig. 2(b)). Line currents from three-phase rectifiers typically contain two positive peaks and two negative peaks ( Fig. 2(c)). Fig. 2(d) and Fig. 2(e) depict the in-phase and quadrature harmonic components contained in the single-phase and three-phase rectified line currents, respectively, as calculated using the Fourier Transform and scaled to maintain the time-domain unit (A). Both contain significant odd-harmonic components, though the three-phase rectified load does not impart 3rd-harmonic components (nor any other zero-sequence harmonic) into the line currents so long as the motor is well balanced. As the speed of the motor varies with the inverted waveform frequency, the current drawn by the motor varies as well. In turn, this causes variations in the harmonic content of the AC line currents, which can thus be used to track the motor's power draw.

IV. THE NILM-CAPABLE ENERGY MONITOR
The NILM-capable energy monitor used to collect data in this study measures the 3-phase voltages and currents supplied to the electrical panel, and a data acquisition system samples the sensor outputs (3 kHz per channel with 16-bit resolution). This high sample rate and resolution allows the precise measurement of both transients useful for identifying individual loads changing states (e.g. turn-on and turn-off) and harmonic content useful in tracking variable loads.
A local computer processes these current and voltage data into their spectral envelopes using the Sinefit algorithm [58], which extracts the spectral envelopes for each phase current over each line cycle. When multiplied by the nominal phase voltage, V ph (e.g. 230 V rms ), these in-phase and quadrature components are, respectively. Here, N represents the number of data points collected over a single line cycle (N is nominally 60 when sampling 50 Hz line-frequency data at a 3kHz sample rate), and k is the line-frequency harmonic (as examples, k = 1 corresponds to 50 Hz, and k = 3 to the 150 Hz component). Of note, P 1 and Q 1 represent the real and reactive power streams under the assumption of constant voltage amplitude and phase. The computer stores these data streams in a NILM-optimized database (NilmDB [59]), from which they can be accessed and manipulated using the NILM Manager and Dashboard software packages [60], [61]. These platforms allow users to create and operate custom load identification algorithms, run equipment and system-level diagnostics, and visualize data, e.g., [41]. As indicated in Figs. 2, the majority of the signal information for variable speed drives is captured in the odd harmonic components below and including the 7th harmonic. As such, we configured the monitor to calculate the real and reactive components for each phase at the fundamental (k = 1), 3rd, 5th, and 7th harmonics (k = 3, 5, 7). From these components, we also configured the monitor to calculate the harmonic components of apparent power,

V. KHALIFA UNIVERSITY RESEARCH FACILITY
The building selected for monitoring is an outlying research facility on the Masdar City Campus of Khalifa University.
The cooling system supplies a total climate controlled area of 228.67 m 2 and features a Toshiba MMY-MAP1204T8 outdoor unit with a rated heat-rejection capacity (q h,rated ) of 33.5 kW. The outdoor unit contains dual rolling-piston compressors and a single condenser fan. Each compressor is electrically powered from a three-phase rectifier/three-phase inverter circuit. The condenser fan is powered from a single-phase rectifier/three-phase inverter circuit. Inside the facility, the refrigeration loop feeds three thermostatically controlled fan-coil evaporator units for direct cooling of zone air and a chilled water heat exchanger for radiant floor cooling. During all tests described in this paper, the refrigerant loop was connected to the chilled water/radiant-floor system.

DIRECT COOLING SYSTEM MONITORING
As part of a previous research project, the cooling system was instrumented to measure, among other things, its rate of heat rejection. Specifically, a thermopile measuring the temperature difference between the condenser inlet and outlet air ( T cond ), and a tachometer measuring the condenser fan speed which was empirically correlated to the volumetric air flow rate across the condenser coils (v a ), provide the measurements required to estimate the outdoor unit's rate of heat rejection,q Here, C p,a is the heat capacity of air and ρ a is the air density, which we estimate as a function of outdoor temperature (T e ) via the ideal gas law, where p a is the absolute air pressure and R a is the specific gas constant for dry air (287.06 J/kg·K).
The instrumentation system also includes temperature sensors on the chilled water system providing chilled water supply and return temperatures (T w,s and T w,r , respectively) and chilled water volumetric flow rate (v w ). These measurements provide estimates of the heat absorbed by the indoor chilled water system as,q where C p,w is the heat capacity of water, ρ w is the water density, and T w = T w,r − T w,s . Finally, current and voltage sensors local to the compressors, condenser fan, and circulating water pump, provide signals from which active power drawn by each of these components, P comps , P cond , and P pump , may be computed. Thus, the installed instrumentation provides two methods for measuring the total heat rejected by the outdoor unit, directly via (6) or indirectly by calculating the total power demand as, and usingq c from (8) to evaluateq h by (2). In general, we estimate the rejected heat as the average of these two calculations.

VI. HEAT REJECTION MAPPING
The heat rejection performance of a properly functioning cooling system can generally be modeled as a nonlinear function of a select few operating variables, e.g., its power draw, the outdoor environmental temperature, and the indoor zonal temperature (T z ) of the air or water at the evaporator [12], [54]. Toshiba provides data relating the heat rejectedq h to the outdoor unit's (compressors plus condenser) power draw, P in,o = P comps + P cond and the outdoor temperature, but only for a single nominal indoor air temperature of 27 • C. Fig. 3 shows these provided data points. The data reveals two characteristics of the cooling system's performance. First, the amount of heat rejected increases approximately logarithmically with input power due to the limited surface area of the heat exchanger. Second, as the outdoor temperature increases, the air conditioner requires more input power to increase the cooling fluid temperature so as to reject a given amount of heat.
This relationship can be empirically modeled as a bi-quadratic function, where, x = ln(P in,o ) and y = T e . The k-terms are constants relating the contribution of the bi-quadratic variables to the heat rejection rate. Fitting this relationship in the least-squares sense to the manufacturer-provided data results in the dashed curves of Fig. 3 and the k-coefficients shown in Table 1.
Prior to the installation of the NILM-capable energy monitor, the performance of the cooling system was extensively tested over a wide range of operating points by varying the speed of the compressors and condenser fan, thus varying the  total power demand of the outdoor unit, and doing so on days and nights spanning a wide range of outdoor temperature. Fig. 4 shows the estimated heat rejection evaluated from data collected during these tests. The dashed lines in these plots are the best-fit curves of Fig. 3. Comparing these color-coded dashed curves with corresponding color-coded data points shows that particularly for times of lower environmental temperatures (the darker and lighter blue data points), the measured heat rejection at a given input power is lower than expected from (10). This is primarily due to low chilled water supply temperatures (T z ) at the evaporator, which decreases the heat rejection performance of the cooling system.
Nominally, the rate of heat transfer in the evaporator is proportional to the difference in temperature between the refrigerant and chilled water. Thus, to account for temperature variations in the chilled water, we first calculate the error in the estimated heat rejection produced by (10) for the data points in Fig. 4, Then, we fit, in the least-squares sense, a linear function of water supply temperature, to , and combine this relationship with the manufacturer defined relationship of (10) to create a heat rejection map that also incorporates the chilled water temperature as an input, i.e., Table 1 also provides the values for the coefficients of (12). Concretely, (13) offsets the manufacturer predicted heat rejection based on the chilled water temperature. Fig. 5 depicts the adjusted data accounting for this offset, showing better matching of the adjusted measured data to the original manufacturer mapping.

VII. ELECTRICAL POWER MAPPING
We installed the NILM-capable energy monitor in the facility's general service subpanel, which distributes electricity to the outdoor unit along with standard building loads including lighting, air handling, and plug loads. Electrical data was collected as described in Sect. IV at the site over a two week period. Fig. 6 shows the outdoor unit's active power (P in,o ) against the total three-phase power measured at the panel by the installed energy monitor. This scatter-plot gives an indication of the outdoor unit's load requirements in relation to the rest of the building. With the unit off, the building load varies between approximately 4 kW and 14 kW. Thus, simply tracking the panel's fundamental power draw does not provide  sufficient information for tracking this variable load. For example, a total panel power draw of 12 kW could correspond to an outdoor unit power draw anywhere between 0 kW and 6 kW.
However, tracking the harmonic content at the panel increases the resolution of the cooling system information by effectively eliminating dominantly linear loads (e.g. fixed-speed motors) and constant power non-linear loads (e.g. lighting ballasts) from the data stream. The bottom figure plots the outdoor unit power against a multilinear regression of select power and harmonic streams measured at the panel. This relationship is far better defined with variability in outdoor unit power for a given fit power now less than 1 kW. Thus, the harmonic streams can be treated as ''features'' from which the cooling system's power draw can be estimated.
Previous studies investigating the relationship between VSD motor loads and harmonic line currents [51], [52] focused on individual harmonic streams as predictors. In this study, we broaden the feature space to include the fundamental real, reactive, apparent, and harmonic streams as defined in (3)-(5) and assess their ability to collectively predict performance in a real-world building environment. To do this, we define a predictor function class,p (H, x), relating a matrix of metered fundamental and harmonic power streams, H, to any directly measured air conditioner power stream, p, through parameter vector, x. These functions can take many forms, e.g. a neural network [62]. In the particular case of the Khalifa University facility's air conditioner system, a multivariate linear model performs remarkably well (in passing, we note that a multivariate linear model is equivalent to a neural network containing nodes with linear activation functions). Thus,p (H, x) = Hx, (14) and the parameter vector can be solved in the least-squares sense as,x With three phases, four harmonics (k = 1, 3, 5, 7) per phase, and three streams (P, Q, S) per harmonic, a matrix H containing all 36 of these features could show multi-collinearity between its columns [63], [64]. Even with large data sets (the analysis of this study utilized approximately 8,000 data points interpolated from the NILM-power streams at 1 min. intervals to match the data collected by the instrumentation system), this multi-collinearity can result in overfitting the data leading to inaccurate coefficient estimates and large prediction errors.
Eliminating unnecessary features from H reduces this problem. One method to do this is to perform the least squares fit across all combinations of predictor variables and use a model quality metric that penalizes the goodness of fit by the number of features used, e.g., adjusted R 2 or Mallow's Cp, to select the best model. However, the 36 harmonic features recorded represent over 68 billion combinations of potential predictor variables making this impractical. Instead, we first structure H based on physical intuition, aggregating some features and eliminating others. and then perform a backwards stepwise regression, iteratively eliminate streams based on confidence estimates in the associated coefficient values. Specifically, we define the feature matrix for the compressors as, 157 ,  a,157 + b,157 + c,157 , 1], (16) and that for the condenser fan as, 157 , 1]. (17) In (16) and (17), the y,1357 term is a matrix with columns containing the panel streams for phase y, i.e., y,1357 = [P y,1 , Q y,1 , S y,1 , P y,3 , Q y,3 , S y,3 , P y,5 , Q y,5 , S y,5 , P y,7 , Q y,7 , S y,7 ]. (18) y,157 is the same matrix but with the 3rd harmonic streams removed, and 1 is a column vector with each element equal to one.
These compressor and condenser fan feature matrices are defined based on the loads' electrical configurations. VOLUME 8, 2020 The compressor is a three-phase load and so should not draw significant three-phase harmonic content as shown in Fig. 2(e). Thus, H comp does not contain any three-phase harmonic features. The condenser fan is a single-phase load powered from phase a. Single-phase non-linear loads do draw third-harmonic currents (Fig. 2(d)) and so H cond contains the third-harmonic content from phase a. The fundamental, 5th, and 7th, harmonic content from phases b and c are included to allow a fit to offset the contributions of three-phase loads (e.g. the compressors) on their corresponding counterparts of phase a.
'Both matrices contain columns aggregating like-harmonic features in a way that maintains the full column space of a matrix of individual features. This is done in preparation for the dimensionality reduction algorithm described below, which we perform to further protect against overfitting the data. Aggregating the like-harmonics effectively averages the individual features and can reduce ''noise'' from extraneous building loads. If however, one harmonic feature correlates better to the power draw of the component of interest, aggregating corrupts the more correlated feature. Constructing these aggregate features while maintaining the column space allows the dimensionality reduction algorithm to select the aggregate or individual columns with parameter coefficients (elements ofx) with the highest relative confidence metric. The individual feature coefficients can then be calculated from the resulting parameter vector estimate,x.

DIMENSIONALITY REDUCTION
The dimensionality reduction algorithm is an iterative approach that eliminates features (columns of H) based on confidence estimates in the associated coefficient values. In this approach, we first set p as P comps or P cond and H to H comps or H cond , respectively. Then, we solve for x as defined in (15), and estimate the residual variance for the fit as, Here, N and M correspond to the dimensions of H (an N × M matrix) and indicate the number of data points observed and feature streams considered, respectively. e = p −p is the fit's residual and so e T e is its sum of squared errors. Withσ 2 , we then generate a variance vector for the estimated coefficients from the diagonal of the coefficient covariance matrix, i.e., Taking the element-wise square root ofσ 2 x gives a standard deviation estimate vectorσ x , whose elements are inversely related to our confidence in their corresponding coefficients ofx. Thus, normalizing this vector by the absolute values of the corresponding coefficients, where indicates element-wise division, provides a confidence metric for each coefficient (this metric is equivalent to the inverse of the t-statistic). The maximum-valued element ofσ x,norm corresponds to the coefficient in which we have the least confidence. We then eliminate the corresponding column in H, effectively eliminating that stream as a feature for determining p, and repeat the process. A variety of criteria can be used to determine when to end this iterative process, e.g., performing a t-test using the inverse of (21). The t-test however assumes normality in the residual distribution, which is not guaranteed here due to the operating schedules of extraneous loads and unmodeled nonlinearities in (14). This can lead to overconfidence in the t-test results. To err on the side of caution, we instead take a graphical, heuristic approach based on the root means square of the error at each iteration,  We perform the iterative process in full, removing columns of H until only the 1 column remains. At each iteration, we evaluate the RMSE, and consider its trend. Fig. 8 and Fig. 9 depict these trends for the fits to the total compressor power draw, P comps , and condenser fan power draw, P cond , respectively. During the beginning of the iterative process, the RMSE does not increase very much indicating that little information is lost by eliminating the initial features. As the number of features decreases towards zero, the RMSE begins to increase more rapidly between iterations, producing a ''knee'' region in the curve where a good tradeoff between fit accuracy and model complexity is often found. We select the illustrated points at the bottom of these knee regions, four features and seven features in the dimension-reduced matrices, H comps,r and H cond,r , respectively.
The corresponding compressor and condenser fan power estimates are defined as, P comps = H comps,rxcomps,r wherex comps,r andx cond,r , are the reduced-dimension coefficient vectors. These estimates can also be defined in terms of only individual harmonic features as, In these equations, the H i matrices contain all of the individual features of the corresponding H r matrices in (23) and (24). That is, H r = H i T, where T is a transformation matrix relating the individual feature columns of H i to the individual and aggregate feature columns of H r . This same transformation matrix, T, relates the individual coefficient vectors of (25) and (26) to the reduced-dimension coefficient vectors of (23) and (24), as,x i = Tx r .   Table 2 shows the reduced-dimension individual features of H comps,i along with the corresponding coefficients of x comps,i , while Table 3 shows those for H cond,i andx cond,i . These coefficient values correspond to the average values obtained during a 10-fold cross validation procedure, where the dimensionality reduction method was performed during each fold. Both H comps,i and H cond,i contained consistent features across all folds. The final columns of these tables give the standard deviations in the estimated coefficient values across the 10 folds. The bottom rows reveal the average RMSE and the standard deviation in RMSE calculated from the model performance. These results show good stability in the estimated coefficient values and stable predictive performance in the both model outputs.    Fig. 10 -Fig. 12 compare the measured power streams of the entire outdoor unit, only the compressor, and only the condenser fan with their corresponding model estimates during VOLUME 8, 2020 a representative three-day testing window. Here, we calculate the model estimate of the outdoor unit power as, The extent to which the estimated power matches the measured power in each plot indicates the ability of the reduced-order harmonic mapping to track the cooling system's component operations. These plots also confirm the dominance of the compressors' power draw (P comps , middle plot) in the outdoor unit's total power demand (P in,o , top plot). However, because the condenser fan is a single-phase load and draws 3rd harmonic components, its smaller power draw (about one-tenth that of the compressors) can still be disaggregated to a reasonable degree from the total power draw at the panel. Of note, during testing the cooling system's operation was controlled manually and not dictated by a thermostat or by any predefined building schedule.

VIII. COMBINED MAPPING PERFORMANCE
Combining the mapping procedures of the previous two sections allows us to estimate the performance of the building's cooling system from the NILM-capable energy monitor. Fig. 13 compares the measured heat rejection with estimates from two combined mappings over a representative 36-hour period. The first estimate, 'Manufacturer Estimate', combines the harmonic map with the heat rejection map of (10) derived solely from manufacturer data, and the second, 'Adjusted Estimate', combines the harmonic map with the adjusted heat rejection map of (13).
At the beginning of this time interval, both estimates match very well with the measured heat rejection, but over time the estimate using only the manufacturer-derived map increasingly overestimates the heat rejection rate. This is due to progressively lower zone and mass temperatures as we tested under the full range of system operating conditions. In turn, the chilled water temperature (Fig. 14), which is not represented in the manufacturer-derived map, decreased significantly over the 36-hour time period. This test represents a ''worst-case'' scenario as the cooling system provides more cooling than normally demanded by the facility. Even so, given measurements of chilled water temperature, the adjusted estimate continues to perform well throughout the test. Thus, supplementing the energy monitor data with outdoor temperature (T e ) and zonal temperature (T z ), which in this case is the chilled water temperature but in direct expansion air conditioners is the indoor air temperature, allows good estimation of a properly functioning cooling system's heat rejection performance over a wide operation range.
While both the outdoor air temperature and the zonal water (or air) temperature could be measured directly and supplied to the NILM-capable energy monitor for the combined performance mapping, this does require additional temperature sensors and thus violates the NILM ''philosophy''. This violation however can be avoided in many situations. First, for most building locations, weather forecasting services can provide accurate local temperature information eliminating the need for direct outdoor temperature measurements. Further, chilled water-based cooling systems like the one used in this study already monitor water temperature, and for direct expansion air conditioners, the zonal temperature is the indoor air temperature, which the thermostat uses for control. As such, a smart thermostat could be configured to make the zonal temperature available to the NILM-capable energy monitor without any further intrusion.

CONDITION MONITORING APPLICATIONS
We have used the combined mapping technique to describe the performance of a properly functioning cooling system. While the identified model won't provide accurate heat rejection estimates during times of faulty operation, it does provide information useful for detecting such faulty operation if combined with direct measurements, or more reasonably, thermally based estimates of the system's heat rejection. As an example, Fig. 15 compares the heat rejection estimates with the directly measured heat rejection during a time when the cooling system operated with low refrigerant levels caused by a previously undetected leak in the refrigerant loop. During this period, these estimates significantly exceeded the measured heat rejection. Thus, comparingq h,map with a thermal model's estimateq h,therm , provides a mechanism for detecting such faulty operation [65]. Here,q h,therm is the output of a thermal model for the facility relating its construction, e.g., building architecture and materials, with environmental conditions such as T e and solar insolation, and the building's operating conditions, e.g., T z and total power use, to estimate the facility's heat rejection profile.

IX. DISCUSSION AND CONCLUSION
The mapping and monitoring techniques presented in this paper exploit the correlation of line harmonics in a building's power demand with the operating performance of its cooling system. By monitoring these line harmonics on each phase, cooling system components with different front-end rectification schemes, i.e. 3-phase rectifiers vs. 1-phase rectifiers, can be disaggregated from each other. The analysis in this paper assumes a linear mapping from these harmonics to the active power draw of the cooling system components. This approach is appropriate for moderately sized facilities where the VSDs of the cooling system are often the dominant nonlinear variable loads. For facilities similar in size to the one studied in this paper, the presence of other VSD loads with a combined power demand of similar or greater magnitude will obfuscate the harmonic-performance relationships of the cooling system. However, it may be possible that with more complex nonlinear mapping techniques, e.g. neural network frameworks [62], multiple large VSD loads can be disaggregated. By enabling single-point monitoring of aggregate variable cooling equipment load, a smart meter or energy monitor that tracks harmonics can help detect equipment faults and even local microclimate problems. LES NORFORD received the B.S. degree in engineering science from Cornell University, in 1973, and the Ph.D. degree in mechanical and aerospace engineering from Princeton University, in 1984. He is currently a Professor with the Department of Architecture, Massachusetts Institute of Technology. He specializes in energy studies, controls, and ventilation and is seeking to improve the way buildings use the earth's resources.