#!/usr/bin/env python# coding: utf-8# This script demonstrates how to use data downloaded from the Intelligent Building Agents Project database located at:# https://www.nist.gov/el/energy-and-environment-division-73200/intelligent-buildings-agents-laboratory-ibal/data/database# # This example focuses on using the data from the lab to examine the performance of the two chillers, Chiller1 and Chiller2.# Both chillers use variable speed scroll compressors. Chiller1 has a nominal capacity of 26.4 kW or 7.5 tons and Chiller2 # has a nominal capacity of 52.8 kW or 15 tons. These chillers replaced the original chillers in October of 2019, # so data before then are excluded from the analysis.# import pandas as pdimport numpy as npimport datetimeimport matplotlib.pyplot as pltdf = pd.read_csv('scaled_data.csv')md = pd.read_csv('MetaData.csv')# Set the index in the metadata as the measurement IDmd.set_index('measurement_id',inplace = True)# Convert the time string to datetime and remove the date column, which is already captured in the time columndf.time = [x[:-6] for x in df.time]df.drop(columns = ['date'],inplace = True)# Format the time string to datetime. Most of the time strings include milliseconds, but some do not. If the milliseconds# are included, the format includes %f at the end.df.time = [datetime.datetime.strptime(x,'%Y-%m-%d %H:%M:%S.%f') if len(x)>19 else \           datetime.datetime.strptime(x,'%Y-%m-%d %H:%M:%S') for x in df.time] # Filter the data to include only times when a chiller is in usedf = df[(df.ch1_on > 0)|(df.ch2_on > 0)]df.index = np.arange(0,np.shape(df)[0])# Convert the chiller setpoints from A to F using the metadatadf.ch1_t_sp = (df.ch1_t_sp - md.loc['ch1_t_sp'].scaling_coeff_a0)/md.loc['ch1_t_sp'].scaling_coeff_a1df.ch2_t_sp = (df.ch2_t_sp - md.loc['ch2_t_sp'].scaling_coeff_a0)/md.loc['ch2_t_sp'].scaling_coeff_a1# Calculate the load on the chiller. This will require the specific heat and density of 30 % propylene glycol.def convertT(T,units = 'IP'):    '''Convert temperature from F to C if units = IP or C to F if units = SI.'''    if (units == 'IP'):        tC = 5/9*(T - 32)        tF = T    else:        tF = T*9/5+32        tC = T    tK = tC + 273.15    return tF,tC,tK # F, C, K    def convertAV(AV,units = 'IP'):    '''Convert the flow rate from gpm to m3/s'''    if units == 'IP':        AV = 0.0000631*AV    else:        AV = AV/0.0000631    return AV # m3/s        def fluidProps(T,units = 'IP',fluid = 'water',conc = 0.3):    '''Calculate properties in liquid. Returns rho, cp, h. conc is concentration with a value between 0 and 1 (1 = 100 %).'''    tF,tC,tK = convertT(T,units)    if (fluid == 'water'):        rho = 995.044+0.278828*tF-0.0046581*tF**2+0.0000219711*tF**3-5.25842*10**(-8)*tF**4        cp = 4.18        h = 0.076174 + 4.18766*tC # [kJ/kg] enthalpy of water at atmospheric pressure; from a curve fit in EES    else: # (fluid == 'pg'):        rho = 508.41109-182.40820*conc+965.76507*273.15/tK+280.29104*conc*273.15/tK-472.22510*(273.15/tK)**2        cp = 4.47642+0.60863*conc-0.71497*273.15/tK-1.93855*conc*273.15/tK+0.47873*(273.15/tK)**2	        h = -20.89 + 3.9095*tC # [kJ/kg] enthalpy of PG at 30 % at atmospheric pressure; from a curve fit in EES    return rho,cp,h # kg/m3, kJ/kg-K, kJ/Kdef mDot(T,AV,fluid,units = 'IP'):    '''Calculate the mass flow rate from the temperature and volumetric flow rate'''    rho,cp,_ = fluidProps(T,'IP','PG')    mDot = rho*convertAV(AV)    return mDot,cp # kg/s, kJ/kg-Kdf['ch1_mDot_e'],df['ch1_cp_e'] = mDot(df.ch1_e_out_rtd,df.ch1_f_e,'PG','IP')df['ch1_mDot_c'],df['ch1_cp_c'] = mDot(df.ch1_c_out_rtd,df.ch1_f_c,'water','IP')df['ch2_mDot_e'],df['ch2_cp_e'] = mDot(df.ch2_e_out_rtd,df.ch2_f_e,'PG','IP')df['ch2_mDot_c'],df['ch2_cp_c'] = mDot(df.ch2_c_out_rtd,df.ch2_f_c,'water','IP')df['q_ch1_e'] = df.ch1_mDot_e*df.ch1_cp_e*(convertT(df.ch1_e_in_rtd)[1]-convertT(df.ch1_e_out_rtd)[1])df['q_ch1_c'] = df.ch1_mDot_c*df.ch1_cp_c*(convertT(df.ch1_c_out_rtd)[1]-convertT(df.pump4_out_rtd)[1])df['q_ch2_e'] = df.ch2_mDot_e*df.ch2_cp_e*(convertT(df.ch1_e_in_rtd)[1]-convertT(df.ch2_e_out_rtd)[1])df['q_ch2_c'] = df.ch2_mDot_c*df.ch2_cp_c*(convertT(df.ch2_c_out_rtd)[1]-convertT(df.pump4_out_rtd)[1])# Calculate the COP for each chillerdf['ch1_cop'] = [a/(b/1000+c/1000) if b > 1000 else 0 for a,b,c in zip(df.q_ch1_e,df.ch1_power,df.pump1_power)]df['ch2_cop'] = [a/(b/1000+c/1000) if b > 1000 else 0 for a,b,c in zip(df.q_ch2_e,df.ch2_power,df.pump2_power)]# Plot some of the results. The plot of COP vs Load shows that Chiller2 is generally more efficient than Chiller1 for the same load.fig, ax = plt.subplots(2,2,figsize = [9,5],sharex = False)ax[0][0].plot(df.index,df.ch1_cop,'.',df.index,df.ch2_cop,'.')ax[0][0].set_xlabel('Index')ax[0][0].set_ylabel('COP')ax[0][0].legend(['CH1','CH2'])ax[0][1].plot(df.q_ch1_e,df.ch1_cop,'o',df.q_ch2_e,df.ch2_cop,'.')ax[0][1].set_xlabel('Load [kW]')ax[0][1].set_ylabel('COP')ax[0][1].legend(['CH1','CH2'])ax[1][0].plot(df.index,df.ch1_t_sp,'.')ax[1][0].set_xlabel('Index')ax[1][0].set_ylabel('Setpoint [F]')ax[0][0].grid()ax[0][1].grid()ax[1][0].grid()ax[1][1].grid()fig.tight_layout()