<jupyter_start><jupyter_text>Create yearly mill usage profile based on measured usageThis notebook will follow a very similar approach to the welder notebook here: https://github.com/RadiantLabs/microgrid-appliance-analysis/blob/master/scripts/usage_profiles/welder_usage_profile.ipynb I will include a lot less context in this notebook for now. I can fill it out later if needed. Definitions**Usage Profile**: A usage profile is a unitless series of how often the mill is used on an hourly basis throughout a year. Being unitless allows us to apply costs, power or other factors to generate load, throughput or cost profiles. **day_hour**: Monday at 2pm is a single `day_hour` for example. Monday at 3pm is a different `day_hour`. There are 24 hours x 7 days = 168 `day_hour`s. General Steps1. Import measured usage data from a mill. This data is in 2 minute increments and the number of mill revolutions every 2 minutes.2. Create a unitless utilization factor based on RPM. Create a separate factor for both kW and grain (production) throughput.3. Resample 2 minute data to 1 hour data, summing new utilization factors. 4. Group data by `day_hour`s and export for web app5. Generate a yearly profile by randomly sampling measured `day_hour`s and applying them to the rest of the year.6. Compare stats of yearly usage profile to measured usage profile to make sure we made reasonable assumptionsWe probably don't need to interpolate, since we have 5-6 measured values per `day_hour` Questions and TODOs* Set column header for datetime (app needs this)* Understand why operators run mills consistently at half RPMs or what the full rated load is:> Amanda: is the engine running harder than full rated load during these times? I think so, because it would then mean that most of the activity is happening around 80% of full rated load, which is normal for a diesel engine* Find performance curves:    1. rpm-to-power: From common mills (https://www.alibaba.com/product-detail/Changfa-type-20-hp-S1115-single_60629667381.html?spm=a2700.7724857.normalList.30.48f321e0kzd6H9), we can see that 16.2kW = 2200 RPM which would imply this mill is running at half speed most of the time. Also, we don't know the kW at 1100 where these mills are mostly running     2. rpm-to-grain-throughput Running this notebookJupyter notebooks are very flexible but have problems with out-of-order function executions.For example, dropping a column once and then trying to drop the same column again becauseyou're evaluating the same cell will throw an error. Go to Cell > Run All from the menu often or if having problems. For final output, always run it.  Setup & Library Imports<jupyter_code># Reset all variables so that we can 'run all cells' and not get unused variables hanging around
%reset -f

# Most of this comes with anaconda distribution. I thinkt the only one 
# you have to install is:
# conda install pivottablejs
import pandas as pd
from scipy.stats import mode
from pivottablejs import pivot_ui
from collections import defaultdict
from functools import partial
import random
import json
import pytz

# Should have pandas 0.23 or greater. If not and you're using Anaconda for packages, 
# do this in the terminal: `conda update pandas`
pd.__version__<jupyter_output><empty_output><jupyter_text>Import Mill Data<jupyter_code># Uncomment / comment based on the mill you want to analyze

## Rice Mill (Original)
# excel_file_path = 'inputs/RiceHulling_LP_Tanzania 20180908-20181018.xlsx'
# sheet_name = 'RiceHulling_LP_Tanzania'
# output_generator_kw_data_name = 'output/rice_mill_kw_generator_data.json'
# output_generator_grain_data_name = 'output/rice_mill_grain_generator_data.json'
# output_profile_name = 'output/rice_mill_usage_profile.csv'

## Rice Mill (UTC)
excel_file_path = 'inputs/RiceHulling_LP_Tanzania 20180908-20181110_UTC.xlsx'
sheet_name = 'download (1)'
output_generator_kw_data_name = 'output/rice_mill_kw_generator_data.json'
output_generator_grain_data_name = 'output/rice_mill_grain_generator_data.json'
output_profile_name = 'output/rice_mill_usage_profile.csv'
# Import Excel file, specify the sheet & import it into a Pandas Dataframe.
# df is short for Pandas DataFrame - it makes it clearer what this datastructure is
df_measured_2min = pd.read_excel(excel_file_path, sheet_name=sheet_name)

# Rename columns so they are shorter and easier to work with.
df_measured_2min = df_measured_2min.rename(columns={'row': 'time_utc', 'value': '2min_count'})

# Keep only relevant columns
df_measured_2min = df_measured_2min[['time_utc', '2min_count']]
df_measured_2min.head()
# (rows, columns)
df_measured_2min.shape<jupyter_output><empty_output><jupyter_text>Convert Timezone<jupyter_code># To see all timezones available (but select only the first 55 to see Africa): 
# pytz.all_timezones[0:55]  
# There was no Tanzania listed. Nairobi is +3 which is the same as Tanzania
# There should not be a problem with daylight savings time - from my research neither observe it
tanzania_tz = pytz.timezone('Africa/Nairobi')
# Convert date string to proper datetime so we can work with timezones
df_measured_2min['time_utc'] = pd.to_datetime(df_measured_2min['time_utc'])
# Add local time
df_measured_2min['time_local'] = df_measured_2min['time_utc'].dt.tz_localize('utc').dt.tz_convert(tanzania_tz)
df_measured_2min = df_measured_2min
df_measured_2min.head()
# Helper functions for making and matching day_hour columns

def shorten_day_name(day_string):
    """Shorten a day name to the first 4 letters (1Saturday => 1sat)
    This requires a string passed in.
    """
    return day_string[0:4].lower()

def composite_val(day_name, hour):
    """Generate a composite string value that can be used for dictionary 
    keys or other uses.
    For example, 1Saturday at 10am => 1sat_10
    """
    padded_hour = str(hour).zfill(2)
    return "{}_{}".format(shorten_day_name(day_name), padded_hour)<jupyter_output><empty_output><jupyter_text>Calculate RPMThe mill sensor measures rotations every 2 minutes. Calculate RPMs by dividing by 2<jupyter_code>def count_to_rpm(count):
    return round(count / 2)

df_measured_2min['rpm'] = df_measured_2min['2min_count'].apply(count_to_rpm)
df_measured_2min.sample(10)<jupyter_output><empty_output><jupyter_text>Explore Data<jupyter_code># The data is skewed because most of it is zeros:
df_measured_2min.describe()
# Filter out all zeros and very low RPMs and look at it again:
# df_measured_nozeros = df_measured_2min[df_measured_2min['rpm'] > -100]  # temp: effectively doesn't filter
df_measured_nozeros = df_measured_2min[df_measured_2min['rpm'] > 100]
df_measured_nozeros.sample(10)
# pivot_ui(df_measured_nozeros, 
#          rows=['2min_count'],
#          rendererName="Table",
#          aggregatorName="Count")
# Plot histogram of rpm including zeros
df_measured_2min.plot.hist(y='rpm', figsize=(20, 4))
# Plot histogram of rpm excluding zeros and low RPMs
df_measured_nozeros.plot.hist(y='rpm', figsize=(20, 4))
# Most common value (mode):
df_measured_nozeros['rpm'].value_counts().idxmax()
# Add day, hour_of_day, day_hour to 2 minute intervals to see the variation in RPMs with different time intervals
# First set index of dataframe to time_local
df_measured_nozeros = df_measured_nozeros.set_index('time_local')
df_measured_nozeros['day'] = df_measured_nozeros.index.dayofweek.map(str) + df_measured_nozeros.index.day_name()
df_measured_nozeros['day'] = df_measured_nozeros["day"].apply(shorten_day_name)
df_measured_nozeros['hour_of_day'] = df_measured_nozeros.index.hour
df_measured_nozeros['day_hour'] = df_measured_nozeros.apply(lambda row: composite_val(row['day'], row['hour_of_day']), axis=1)
df_measured_nozeros['week'] = df_measured_nozeros.index.week
df_measured_nozeros['rpm_binned'] = pd.cut(df_measured_nozeros['rpm'], 5, precision=0, )
df_measured_nozeros.sample(10)
pivot_ui(df_measured_nozeros, 
         rows=['rpm_binned'],
         cols=['hour_of_day'],
         rendererName="Bar Chart",
         aggregatorName="Count")
pivot_ui(df_measured_nozeros, 
         rows=['rpm_binned'],
         cols=['day'],
         rendererName="Bar Chart",
         aggregatorName="Count")<jupyter_output><empty_output><jupyter_text>Thoughts* We're seeing the RPM average and mode as ~1000. We think the max RPM is 2200. We are measuring rotations every 2 minutes but 2 minutes is fairly long. Is it possible these mills cycle on and off (or low rpms) within those 2 minutes? For example, they could run at max for 30 seconds, spin down to low RPMs for 30, then spin back up to max continuously?     1. The data could support that type of cycling if that behavior was programmed into the mill    2. The data is too consistent to support people manually cycling these mills. I would expect more variance in the distributions of RPMs.  Approach to convert RPM to unitless usage countsRPMs relate to 3 quantities that we care about:* kW* kWh* Grain throughputThese quantities are directly related to costs and revenue for the appliance owner and grid operator. We have measured revolution count in 2-minute intervals. The average RPM is half the number of mill revolutions.However, if we want to work with hourly intervals, then we can't add up RPMs to get the total RPMs within an hour. We have a few options:1. Average the RPMs across the hour. If the mill runs for 30 minutes at 1200 RPMs, then the average would be 600 RPMs. Grain throughput and kW likely have a non-linear relationship to RPMs, so this would really distort the results.2. Assign a unitless usage factor to each minute interval. If the mill is running at full production (throughput and kW), then we can assign it a value of 1 (2 for a 2-minute interval). We can scale the factor down as the RPMs go down. The second approach is what I will use. Throughput vs. kWTo calculate kW (and therefore kWh) and throughput, we need to know how RPM relates to these two factors. There is no reason to assume the kW and grain throughput vs RPM scale the same. So let's create two functions that calculate kW and grain throughput independently.For now, these factors will be the same. But the functions can be easily modified as we get more information without refactoring code downstream.  RPM to unitless utilization factorsThis mill is likely rated at 16.2kW @ 2200 RPM. Let's assume some function that relates it's utilization given an RPM.<jupyter_code># What we have measured from the mill sensors is revolutions per 2 minutes: `rp2m`. (At this point we've already converted to rpm by dividing it by 2)
# What we want is the percent that this mill is fully utilized over the course of an hour: `kw_factor`.
# This will allow us to multiply this factor by the appliance power (nominal power in kW) to get the kW for that hour.
# This should stay true as long as kW scales linearly with RPM. We can refine later as we get more data.

# We may want to filter very low RPMs out. A slow motor still uses power and mills grain, but I don't know the cutoff
# Assuming we have already 
def rpm_to_kw_utilization(rpm, full_capacity=1110):
    """
    Convert RPM to a kW utilization factor. 
    This assumes we have already converted the 2-minute revolution count (rp2m) count to 1 minute (rpm) 
    """
    # Fraction of full capacity over 2 minutes 
    two_min_rpm_full_capacity_percent = rpm / full_capacity
    
    # There are 30 2-minute intervals in an hour
    two_min_intervals_per_hour_count = 30
    
    # Calculate the fraction of full utilization for the *entire hour* this 2 minute interval
    # contributes. Then when we resample fromm 2min to 1 hour intervals, these will sum correctly.
    # This is a little unintuitive - we could divide by 30 after resampling, but then this 
    # calculation is spread out across this notebook and may be more confusing
    return round(two_min_rpm_full_capacity_percent / two_min_intervals_per_hour_count, 5)
    
rpm_to_kw_utilization(1121)
# Until we have better data, assume the utilization factors are the same between kW and grain
# See explanation above

def rpm_to_grain_utilization(rpm, full_capacity=1110):
    """
    Convert RPM to a grain milled utilization factor
    """
    # Fraction of full capacity over 2 minutes 
    two_min_rpm_full_capacity_percent = rpm / full_capacity
    
    # There are 30 2-minute intervals in an hour
    two_min_intervals_per_hour_count = 30
    
    # Calculate the fraction of full utilization for the *entire hour* this 2 minute interval
    # contributes. Then when we resample fromm 2min to 1 hour intervals, these will sum correctly.
    # This is a little unintuitive - we could divide by 30 after resampling, but then this 
    # calculation is spread out across this notebook and may be more confusing
    return round(two_min_rpm_full_capacity_percent / two_min_intervals_per_hour_count, 5)

rpm_to_grain_utilization(1121)
df_measured_2min['kw_factor'] = df_measured_2min['rpm'].apply(rpm_to_kw_utilization) 
df_measured_2min['production_factor'] = df_measured_2min['rpm'].apply(rpm_to_grain_utilization) 

# Show a single hour: 
df_measured_2min[30220:30250]
# df_measured_2min[30220:30250].to_csv('test_summation', index_label='index')
# Sum the kw_factor for this hour: 2018-10-19 16:00 (axis=0 is column)
# You can see the utilization kw_factor is 0.87098. This is the same as that same hour after resampling to hourly
df_measured_2min[30220:30250].sum(axis=0)<jupyter_output><empty_output><jupyter_text>Resample 2 minute to hourly intervals<jupyter_code># First we need to set the dataframe index to be the time_local datetime
df_measured_2min_index = df_measured_2min.set_index('time_local')
df_measured_2min_index.shape
# Resample while summing every 2-min interval within an hour ('H')
# Drop the now-unused 2min_count
# UTC time will automatically be dropped since you can't sum it
df_measured = df_measured_2min_index.resample('H').sum().drop(columns=['2min_count'])
df_measured.shape  # (rows, columns) where column count doesn't include the index
# This should match the df_measured_2min[30220:30250].sum(axis=0) value for kw_factor above
df_measured['2018-10-19 16:00':'2018-10-19 16:00']<jupyter_output><empty_output><jupyter_text>Add hour, day of week, day_hour columns These columns will be used later for generating yearly usage profile and aggregate stats.<jupyter_code># Add the name of the day of the week to the dataframe (Saturday).
# Prepend that name with a number of the day of the week.
# Monday is 0, Tuesday is 1 and so on. This will allows tools to 
# order the days so they are in order: 0Monday, 1Tuesday, otherwise 
# they will be ordered alphabetical.
df_measured["day"] = df_measured.index.dayofweek.map(str) + df_measured.index.day_name()
df_measured["day"] = df_measured["day"].apply(shorten_day_name)

# Add hour of day (as a number)
df_measured['hour_of_day'] = df_measured.index.hour 

# Add day_hour. For example: 4fri_10
# Possible source of confusion: 
# 4fri is just friday. 4fri_10 is Friday at 10am. 
df_measured["day_hour"] = df_measured.apply(lambda row: composite_val(row['day'], row['hour_of_day']), axis=1)

df_measured.sample(10)<jupyter_output><empty_output><jupyter_text>Generate data for yearly usage profile sampling<jupyter_code>def create_usage_profile_data(df, column):
    """
    Create a dictionary, where each key is a day_hour and each value 
    is a list of measured usage factor values. 
    Takes a Pandas dataframe and returns a python dictionary that can be 
    encoded into JSON for other applications.
    We may be able to use groupby for a more succinct function, but this works
    """
    dict = defaultdict(list)
    for index, row in df.iterrows():
        key = row['day_hour']
        dict[key].append(row[column])
    return dict
measured_usage_kw_factors = create_usage_profile_data(df_measured, 'kw_factor')
# measured_usage_kw_factors
measured_usage_grain_factors = create_usage_profile_data(df_measured, 'production_factor')
# measured_usage_grain_factors<jupyter_output><empty_output><jupyter_text>Export usage data for web app<jupyter_code>with open(output_generator_kw_data_name, 'w') as fp:
    json.dump(measured_usage_kw_factors, fp)
    
with open(output_generator_grain_data_name, 'w') as fp:
    json.dump(measured_usage_grain_factors, fp)<jupyter_output><empty_output><jupyter_text>Generating yearly usage profileNow that we have usage profile data, generate a complete year's usage profile<jupyter_code>def create_year_range_df(year=2018):
    """
    Creates a dataframe with a full year's dates as the index. 
    Add extra derived columns based on that datetime index: 
    (hour, day, hour_of_day, day_hour).
    This dataframe does not contain any appliance data
    """
    start_date_str = '1/1/{}'.format(year + 1)
    start_date = pd.to_datetime(start_date_str) - pd.Timedelta(days=365)
    hourly_periods = 8760
    date_range = pd.date_range(start_date, periods=hourly_periods, freq='H')
    year_hours = list(range(len(date_range)))
    
    # Create a full year with a datetime index (8760 hours)
    df_year = pd.DataFrame({"hour": year_hours}, index=date_range)
    
    # Now add day of week, hour of day and day_hour columns
    df_year['day'] = df_year.index.dayofweek.map(str) + df_year.index.day_name()
    df_year['day'] = df_year["day"].apply(shorten_day_name)
    df_year['hour_of_day'] = df_year.index.hour
    df_year["day_hour"] = df_year.apply(lambda row: composite_val(row['day'], row['hour_of_day']), axis=1)
    return df_year

# Uncomment these to test results. 
# This function is called from generate_usage_profile() 
# df_year_example = create_year_range_df()
# df_year_example.head()
#df_year['kw_factor'] = df_year.apply(partial(sample_usage_index, kw_factor), axis=1)
# df_year['production_factor'] = df_year.apply(partial(sample_usage_index, production_factor), axis=1) 


def sample_usage_index(usage_factors, day_hour):
    """
    Takes the measured usage dictionary and a dataframe row
    from the empty yearly profile created in create_year_range_df. 
    Using the day_hour from that dataframe row, take a random
    sample of the same day_hour from the measured data.
    """
    choice_index_range = range(len(usage_factors[day_hour]))
    return random.choice(choice_index_range)
    

def generate_usage_profile(kw_factors, production_factors, year=2018):
    """
    Generate a full year's usage profile by sampling measured data 
    Takes the measured usage dictionary for both kW and grain factor
    and optional year.
    Returns a dataframe of hourly intevals with sampled appliance values
    """
    # Create dataframe with a full year of hourly intervals
    df_year = create_year_range_df(year)

    # We need to make sure we take the same sample from both kw_factors and production_factors
    # So instead of returning a factor from the sampling, return an index of the list of factors  
    for row_index, row in df_year.iterrows():
        day_hour = row['day_hour']
        usage_index = sample_usage_index(kw_factors, day_hour)
        
        # Add new columns with the sample values
        df_year.at[row_index, 'kw_factor'] = kw_factors[day_hour][usage_index]
        df_year.at[row_index, 'production_factor'] = production_factors[day_hour][usage_index]
        
    return df_year


df_generated_usage_profile = generate_usage_profile(measured_usage_kw_factors, measured_usage_grain_factors)
df_generated_usage_profile.head(10)
df_generated_usage_profile.describe()
# Measured Hourly Averages
pivot_ui(df_measured,
        cols=['hour_of_day'],
        rendererName='Table Barchart',
        aggregatorName='Average',
        vals=['kw_factor'])
# ![Screenshot](screenshots/rice_mill_measured_hourly_avg.png)
# Generated Hourly Averages
pivot_ui(df_generated_usage_profile,
        cols=['hour_of_day'],
        rendererName='Table Barchart',
        aggregatorName='Average',
        vals=['kw_factor'])
# ![Screenshot](screenshots/rice_mill_generated_hourly_avg.png)
# ![Screenshot](screenshots/rice_mill_measured_daily_avg.png)
# Measured Daily Averages
pivot_ui(df_measured,
        cols=['day'],
        rendererName='Table Barchart',
        aggregatorName='Average',
        vals=['kw_factor'])
# Generated Daily Averages
pivot_ui(df_generated_usage_profile,
        cols=['day'],
        rendererName='Table Barchart',
        aggregatorName='Average',
        vals=['kw_factor'])
# ![Screenshot](screenshots/rice_mill_generated_daily_avg.png)
# Show a week's worth of mill usage:
df_week_6 = df_generated_usage_profile.loc[df_generated_usage_profile.index.week == 5]
df_week_6.plot(y='kw_factor', figsize=(20, 4))
# Show a month's worth of mill usage:
df_february = df_generated_usage_profile.loc[df_generated_usage_profile.index.month == 4]
df_february.plot(y='kw_factor', figsize=(20, 4))
# Show entire year's worth of mill usage:
df_generated_usage_profile.plot(y='kw_factor', figsize=(20, 4))<jupyter_output><empty_output><jupyter_text>Export yearly usage profile The web app doesn't need this data but it can be used for other analysis<jupyter_code>print("outputting to: ", output_profile_name)
df_generated_usage_profile.to_csv(output_profile_name, index_label='datetime')<jupyter_output>outputting to:  output/rice_mill_usage_profile.csv
{'fraction_non_alphanumeric': 0.07645979391144779, 'fraction_numerical': 0.019856020326542136, 'mean_word_length': 5.003954802259887, 'pattern_counts': {'":': 1, '<': 48, '<?xml version=': 0, '>': 54, 'https://': 2, 'lorem ipsum': 0, 'www.': 1, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '12', 'dir': 'jupyter-structured-clean-dedup', 'id': '91568', 'license': 'Apache-2.0', 'n_tokens_mistral': 7078, 'n_tokens_neox': 6403, 'n_words': 2274, 'path': 'usage_profiles/rice_mill_usage_profile.ipynb', 'repo_name': 'RadiantLabs/microgrid-appliance-usage-profile-generators', 'size': '21253'}