<jupyter_start><jupyter_text>Sunset solarWe would like to understand how existing supply and demand on the grid changes the marginal value of new solar in the system.This requires us to understand:- The price of electricity at different times of day- The different generation profiles of solar panels at different aspectsIf the current price doesn't support a different aspect of solar panel, we might add- The carbon intensity of the grid at different times of day and thus- The carbon price that would make West facing solar panels more profitable than South facing ones<jupyter_code>import pandas as pd
import plotly.express as px
import numpy as np
from numpy import pi<jupyter_output><empty_output><jupyter_text>I think the scale of this data is GBP/MWH<jupyter_code>price_df = pd.read_csv('DA 2019.csv', parse_dates=['settlement_period']).rename({'N2EX_Index__PRICE': 'Day Ahead Price', 'settlement_period': 'Timestamp'}, axis=1)
px.line(price_df, x='Timestamp', y='Day Ahead Price')<jupyter_output><empty_output><jupyter_text>Remove the outlier<jupyter_code>price_df = price_df[price_df['Day Ahead Price'] < 150]
price_df['Hour'] = price_df['Timestamp'].apply(lambda x: x.hour)
price_df['Month'] = price_df['Timestamp'].apply(lambda x: x.month)

gb = price_df.groupby('Hour').mean()
fig = px.area(gb, x=gb.index, y='Day Ahead Price', title='Daily Variation in Energy Prices: The Double Hump')
fig.update_layout(yaxis=dict(range=[0, 60]))

gb = price_df.groupby('Month').mean()
px.area(gb, x=gb.index, y='Day Ahead Price', title='Monthly Variation in Energy Prices: January is expensive')<jupyter_output><empty_output><jupyter_text>Solar irradiance dataData is provided from the [CEDA Archive](http://data.ceda.ac.uk/badc/ukmo-midas-open/data/uk-radiation-obs/dataset-version-202007/avon/00676_filton/qc-version-1) as horizontal irradiance. Luckily Steph previously wrote some code to convert this according to aspect (referred to as Azimuth in solar talk). We'll use some data from the Bristol area, Filton.<jupyter_code>sun_df = pd.read_csv('midas-open_uk-radiation-obs_dv-202007_avon_00676_filton_qcv-1_2017.csv', header=75, parse_dates=['ob_end_time'])
sun_df = sun_df.rename({'ob_end_time': 'Timestamp', 'glbl_irad_amt': 'Value'}, axis=1)
px.line(sun_df, x='Timestamp', y='Value')<jupyter_output><empty_output><jupyter_text>Looks to me like there are a couple of outliers to be removed. There is a spike on the last day of every month, just before midnight - is this the monthly total or something?<jupyter_code>sun_df = sun_df[sun_df['id_type']=='DCNN']
px.line(sun_df, x='Timestamp', y='Value')<jupyter_output><empty_output><jupyter_text>Much better :)  Data refers to the _previous_ hour. Settlement periods are aligned to the start of the period so we should shift everything by an hour.<jupyter_code>sun_df = sun_df.set_index('Timestamp').shift(-1).reset_index()
sun_df["Timestamp"] = sun_df["Timestamp"].apply(pd.Timestamp)
sun_df['Hour']= sun_df['Timestamp'].apply(lambda x: x.hour)
sun_df['Month'] = sun_df['Timestamp'].apply(lambda x: x.month)

gb = sun_df.groupby('Hour').sum()
px.line(gb, x=gb.index, y='Value', title='Total Horizontal Irradiance by Hour')
gb = sun_df.groupby('Month').sum()
px.line(gb, x=gb.index, y='Value', title='Total Horizontal Irradiance by Month')<jupyter_output><empty_output><jupyter_text>Sun angle and heightIn addition to the solar irradiance on the horizontal, we need to know the locaton (height and azimuth) of the sun throughout the year.Luckily Steph already figured this out and we can get the data here: https://susdesign.com/sunposition/I've downloaded it in _clock time_<jupyter_code>sun_position = pd.read_csv('bristol_sun_angles.csv')
sun_position.drop(0, axis=0, inplace=True)
sun_position.head()
sun_position['Timestamp'] = sun_position.apply(lambda row: pd.Timestamp(f"{row['date']}/2017 {row['clock time']}"), axis=1)
sun_position.drop(['date', 'clock time'], axis=1, inplace=True)

# Nicked from Steph
sun_position[['altitude', 'azimuth']] = sun_position[['altitude', 'azimuth']].apply(np.deg2rad)
px.line(sun_position, x='Timestamp', y=['altitude','azimuth'])<jupyter_output><empty_output><jupyter_text>Fold together irradiance and azimuth data Some functions stolen from Steph<jupyter_code># Resample hourly first
sun_df = sun_position.set_index('Timestamp').resample('H').mean().join(sun_df.set_index('Timestamp')['Value'])
px.line(sun_df, x=sun_df.index, y=['Value','azimuth','altitude'])<jupyter_output><empty_output><jupyter_text>We derive the intensity on a vertical surface at a given azimuth in two steps:1. Turn horiztonal irradiance to vertical irradiance of a plane facing the sun perfectly2. Weight that according to the azimuth  1. Irradiance on a vertical and tilted surface Called "sunflower" intensity, rather charmingly. This describes the fact that it rotates to face the sun. Average UK home has a tilt from 40-50 degrees.<jupyter_code>def horizontal_to_vertical_weighting(altitude, min_alt=np.deg2rad(2)):
    weighting = 1 / np.tan(altitude)
    return weighting if altitude > min_alt else 0

def horizontal_to_inclined(altitude, incline_to_horizontal_deg=45, min_alt=np.deg2rad(2)):
    # Steph's version
#     weighting = np.sin(altitude + incline_to_horizontal_deg/180)/np.sin(altitude)
    
    # Corrected
    weighting = np.sin(altitude + np.deg2rad(incline_to_horizontal_deg))/np.sin(altitude)
    
    # Archy version
#     weighting = (1 / np.tan(altitude)) / np.cos(np.deg2rad(incline_to_horizontal_deg))
    return weighting if altitude > min_alt else 0

sun_df['weighting_vertical'] = sun_df['altitude'].apply(horizontal_to_vertical_weighting)
sun_df['sunflower_vertical'] = sun_df['weighting_vertical'] * sun_df['Value']

sun_df['weighting_tilted'] = sun_df['altitude'].apply(horizontal_to_inclined)
sun_df['sunflower_tilted'] = sun_df['weighting_tilted'] * sun_df['Value']
px.line(sun_df.loc['2017-06-06'], x=sun_df.loc['2017-06-06'].index, y=['weighting_tilted', 'weighting_vertical'],
       title='Vertical surfaces get more sun in the morning and evening, but more horizontally tilted ones benefit at midday')
px.line(sun_df.loc['2017-06-06'], x=sun_df.loc['2017-06-06'].index, y=['sunflower_tilted', 'sunflower_vertical'],
       title='Multiplied through by measured intensity shows benefit of having a more horizontal panel')<jupyter_output><empty_output><jupyter_text>2. Irradiance by azimuth<jupyter_code># Nicked from Steph
orientation_angles = {
    'N': (12/16*(2*np.pi), '#d62728'),
    'NNE': (13/16*(2*np.pi),'rebeccapurple'),
    'NE': (14/16*(2*np.pi),'darkgreen'),
    'ENE': (15/16*(2*np.pi), 'black'),
    'E': (0*(2*np.pi),'#2ca02c'),
    'ESE': (1/16*(2*np.pi),'magenta'),
    'SE': (2/16*(2*np.pi),'midnightblue'),
    'SSE': (3/16*(2*np.pi),'grey'),
    'S': (4/16*(2*np.pi),'#ff7f0e'),
    'SSW': (5/16*(2*np.pi),'dodgerblue'),
    'SW': (6/16*(2*np.pi),'cyan'),
    'WSW': (7/16*(2*np.pi),'darkred'),
    'W': (8/16*(2*np.pi),'#1f77b4'),
    'WNW': (9/16*(2*np.pi),'crimson'),
    'NW': (10/16*(2*np.pi),'lime'),
    'NNW': (11/16*(2*np.pi),'red'),
    }<jupyter_output><empty_output><jupyter_text>For each possible orientation, we need to take angle of incidence of the sun. If the angle is 90 then the sun is directly hitting the surface and we're capturing all of the sunlight (i.e. we have the sunflower intensity). If not, we need to take `sin(x)`where x is the angle of attack<jupyter_code>assert orientation_angles['S'][0] == np.pi/2
def sunflower_to_azimuth(azimuth, plane_orientation):
    """
    Work out the % of sun that falls at a given azimuth
    """
    angle =  plane_orientation - azimuth
    weighting = np.sin(angle)
    if weighting < 0:
        # This corresponds to sun on the other side of the surface
        weighting = 0
    return weighting
sun_df['south weighting'] = sun_df['azimuth'].apply(lambda x: sunflower_to_azimuth(x, orientation_angles['S'][0]))
sun_df['west weighting'] = sun_df['azimuth'].apply(lambda x: sunflower_to_azimuth(x, orientation_angles['W'][0]))
sun_df['north weighting'] = sun_df['azimuth'].apply(lambda x: sunflower_to_azimuth(x, orientation_angles['N'][0]))
sun_df['east weighting'] = sun_df['azimuth'].apply(lambda x: sunflower_to_azimuth(x, orientation_angles['E'][0]))<jupyter_output><empty_output><jupyter_text>Verify that this looks as we expect:<jupyter_code>px.line(sun_df.loc['2017-08-01'], x=sun_df.loc['2017-08-01'].index, y=['south weighting', 'west weighting'],
       title='South facing panels weighted higher at midday, West in the evening')<jupyter_output><empty_output><jupyter_text>Fold them togetherMultiply through the vertical weightings by sunflower<jupyter_code>for direction in ['south', 'west', 'north', 'east']:
    sun_df[f'{direction} intensity'] = sun_df[f'{direction} weighting'] * sun_df['sunflower_tilted']
px.line(sun_df.loc['2017-08'], x=sun_df.loc['2017-08'].index, y=['south intensity', 'west intensity'])
sun_df['Hour'] = sun_df.reset_index()['Timestamp'].apply(lambda x: int(x.hour)).values
sun_df['Month'] = sun_df.reset_index()['Timestamp'].apply(lambda x: int(x.month)).values
gb = sun_df.groupby('Hour').mean().reset_index().sort_values(by='Hour')
px.line(gb, x=gb['Hour'], y=['south intensity', 'west intensity'], title='Averaged over the year')
gb = sun_df.groupby('Month').mean().reset_index().sort_values(by='Month')
px.line(gb, x=gb['Month'], y=['south intensity', 'west intensity'], title='South facing panels have better average monthly performance every month')<jupyter_output><empty_output><jupyter_text>Multiply solar profiles by month * hour averages<jupyter_code>gb = price_df.groupby(['Month', 'Hour']).mean().reset_index()
fig = px.bar(gb, y='Day Ahead Price', x='Hour', facet_row='Month', 
             facet_row_spacing=0.01, height=1000, width=800, 
             color='Day Ahead Price', color_continuous_scale='BuPu')

# hide and lock down axes
fig.update_xaxes(visible=False, fixedrange=True)
fig.update_yaxes(visible=False, fixedrange=True)
# strip down the rest of the plot
fig.update_layout(
    showlegend=False,
    plot_bgcolor="white",
    margin=dict(t=10,l=10,b=10,r=10),
    coloraxis=None,
)



# disable the modebar for such a small plot
fig.show(config=dict(displayModeBar=False))<jupyter_output><empty_output><jupyter_text>Perhaps easier to see as a matrix:<jupyter_code>px.imshow(gb.pivot_table(values='Day Ahead Price', index=['Month'], columns=['Hour']))<jupyter_output><empty_output><jupyter_text>Join average price data<jupyter_code>sun_df = sun_df.reset_index().set_index(['Month', 'Hour'])
avg_prices = price_df.groupby(['Month', 'Hour']).mean()
df = sun_df.join(avg_prices, how='left').sort_values(by='Timestamp')
df
df['south earnings'] = df['south intensity'] * df['Day Ahead Price']
df['west earnings'] = df['west intensity'] * df['Day Ahead Price']
px.line(df, x='Timestamp', y=['south earnings', 'west earnings'], title='This does not look good for Western panels')
gb = df.groupby('Month').sum()
px.bar(gb, x=gb.index, y=['south earnings', 'west earnings'], barmode='group',
      title='A 45 degree S panel beats a 45 W panel every month of the year')
gb = df.groupby('Hour').sum()
px.bar(gb, x=gb.index, y=['south earnings', 'west earnings'], barmode='group',
      title='Even though a W panel does get more evening light throughout the year')
px.scatter(df, x='Day Ahead Price', y=['south intensity', 'west intensity'])<jupyter_output><empty_output><jupyter_text>Finding the optimal We'll do a little grid search over orientation and angles<jupyter_code>def calculate_yearly_income(df, orientation, tilt):
    
    df['tilt_weighting'] = df['altitude'].apply(lambda x: horizontal_to_inclined(x, incline_to_horizontal_deg=tilt))
    df['orientation_weighting'] = df['azimuth'].apply(lambda x: sunflower_to_azimuth(x, orientation))
    
    df['output_intensity'] = df['Value'] * df['tilt_weighting'] * df['orientation_weighting']
    df['cash'] = df['output_intensity'] * df['Day Ahead Price']
    
    return df
    
money = {}
sun = {}
for orientation, (azimuth, _) in orientation_angles.items():
    money[orientation] = {}
    sun[orientation] = {}
    for tilt in range(0, 90, 10):
        output = calculate_yearly_income(df, orientation=azimuth, tilt=tilt)
        
        money[orientation][tilt] =  output["cash"].sum()
        sun[orientation][tilt] =  output["output_intensity"].sum()
        
        print(orientation, tilt)
cash_df = pd.DataFrame(money)
sun_df = pd.DataFrame(sun)
px.imshow(sun_df, title='South facing 50 gets the most sun')
px.imshow(cash_df, title='SSE 60 degrees are the most profitable')<jupyter_output><empty_output><jupyter_text>So this is not a good idea- The loss of energy at other times of day does not justify the increased energy you get in the evening spots- The morning hump is much better served by the South facing panel!  What about carbon intensity?The reason the evening prices aren't higher is that the grid become carbon intense with winter peakers. What if we had a carbon price?Data from https://data.nationalgrideso.com/carbon-intensity1/historic-generation-mix<jupyter_code>grid_intensity = pd.read_csv('df_fuel_ckan.csv', parse_dates=['DATETIME']).set_index("DATETIME")
grid_intensity = grid_intensity.loc['2019']
grid_intensity['Hour'] = grid_intensity.reset_index()['DATETIME'].apply(lambda x: pd.Timestamp(x).hour).values
grid_intensity['Month'] = grid_intensity.reset_index()['DATETIME'].apply(lambda x: pd.Timestamp(x).month).values

px.line(grid_intensity,grid_intensity.index, 'CARBON_INTENSITY')
gb = grid_intensity.groupby('Hour').mean()
px.line(gb, x=gb.index, y='CARBON_INTENSITY', title='CO2 intensity mimics price')
gb = grid_intensity.groupby('Month').mean()
px.line(gb, x=gb.index, y='CARBON_INTENSITY', title='CO2 intensity mimics price')<jupyter_output><empty_output><jupyter_text>Grid intensity is g carbon per kwh.  This report suggests a price of 75 GBP / tonne in 2030 might be reasonable: https://www.edie.net/news/6/EU-carbon-price-set-to-rise-to-EUR32-by-2030--but-experts-say-EUR81-necessary-to-achieve-net-zero-in-the-UK/<jupyter_code>price_per_ton_co2 = 75
price_per_gram_co2 = price_per_ton_co2/1e6
kwh_to_mwh = 1/1e3
grid_intensity['Carbon Price GBP/MWh'] = grid_intensity['CARBON_INTENSITY']*price_per_gram_co2 / kwh_to_mwh<jupyter_output><empty_output><jupyter_text>Merge our grid intensity onto our prices<jupyter_code>carbon_price_df = price_df.set_index('Timestamp').join(grid_intensity[['CARBON_INTENSITY', 'Carbon Price GBP/MWh']])
carbon_price_df['Hour'] = carbon_price_df.reset_index()['Timestamp'].apply(lambda x: pd.Timestamp(x).hour).values
carbon_price_df['Month'] = carbon_price_df.reset_index()['Timestamp'].apply(lambda x: pd.Timestamp(x).month).values
px.line(carbon_price_df, x=carbon_price_df.index, y=['Day Ahead Price', 'Carbon Price GBP/MWh'],
       title='A carbon price would be a substantial portion of cost')
gb = carbon_price_df.groupby('Hour').sum()
px.bar(gb, x=gb.index, y=['Day Ahead Price', 'Carbon Price GBP/MWh'], barmode='group')<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.14675609270203865, 'fraction_numerical': 0.019390397768776147, 'mean_word_length': 6.328467153284672, 'pattern_counts': {'":': 0, '<': 83, '<?xml version=': 0, '>': 83, 'https://': 3, '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': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '436607', 'license': 'MIT', 'n_tokens_mistral': 5407, 'n_tokens_neox': 5088, 'n_words': 1217, 'path': 'sunset_solar.ipynb', 'repo_name': 'archydeberker/sunset-solar', 'size': '15059'}