%load_ext watermark
%watermark -a 'cs224' -u -d -v -p numpy,scipy,pandas,matplotlib,pvlib
Author: cs224 Last updated: 2023-08-16 Python implementation: CPython Python version : 3.11.4 IPython version : 8.12.0 numpy : 1.24.3 scipy : 1.10.1 pandas : 1.5.3 matplotlib: 3.7.1 pvlib : 0.10.1
%matplotlib inline
import numpy as np, scipy, pandas as pd, matplotlib.pyplot as plt, seaborn as sns
import matplotlib.patches, matplotlib.transforms, matplotlib.dates
import tzlocal
import pvlib, pvlib.iotools, pvlib.modelchain, pvlib.location, pvlib.pvsystem, pvlib.temperature, pvlib.irradiance, pvlib.solarposition
pd.set_option('display.max_columns', 500)
pd.set_option('display.width', 1000)
# pd.set_option('display.float_format', lambda x: '%.2f' % x)
np.set_printoptions(edgeitems=10)
np.set_printoptions(linewidth=1000)
np.set_printoptions(suppress=True)
np.core.arrayprint._line_width = 180
SEED = 42
np.random.seed(SEED)
sns.set()
from IPython.display import display, HTML
display(HTML("<style>.container { width:70% !important; }</style>"))
This jupyter notebook belongs to the blog post photovoltaics-balcony-power-station.
The below values will differ from the values in the blog post, because I changed the latitude and longitude.
latitude = 48.1387900033257
longitude = 11.573658530854429
altitude=350.0
balcony_length=3.547
balcony_width=2.031
module_length=1.692
module_width=1.029
angle = 23.0
delta_x = 0.02
delta_y = 0.10
start = pd.Timestamp('2023-08-11 12:30').tz_localize(tzlocal.get_localzone()) # until 16:30
time_slice_minutes = 5
periods = 4*(60/time_slice_minutes)+1
freq=f'{time_slice_minutes}T'
%matplotlib inline
plt.figure(figsize=(15, 8), dpi=80, facecolor='w', edgecolor='k')
ax = plt.subplot(1, 1, 1)
ax.set_aspect('equal')
transform = matplotlib.transforms.Affine2D().rotate_deg(angle)
trans_data = transform + ax.transData
balkony = matplotlib.patches.Rectangle((0.0, 0.0), balcony_length, balcony_width, alpha=0.1, facecolor='black')
balkony.set_transform(trans_data)
ax.add_patch(balkony)
pv_module_1 = matplotlib.patches.Rectangle((0.0, 0.0), module_width, module_length, alpha=0.2, facecolor='yellow')
pv_module_1.set_transform(trans_data)
ax.add_patch(pv_module_1)
pv_module_2 = matplotlib.patches.Rectangle((module_width, 0.0), module_width, module_length, alpha=0.2, facecolor='green')
pv_module_2.set_transform(trans_data)
ax.add_patch(pv_module_2)
origin = np.array([0, 0])
yhat = np.array([0, 1])
x, y, arrow_length = 0.5, 0.5, 0.1
ax.annotate('N', xy=(x, y), xytext=(x, y-arrow_length), arrowprops=dict(facecolor='black', width=5, headwidth=15), ha='center', va='center', fontsize=20, xycoords=ax.transAxes)
ax.set_xlim([-1.2, 3.3])
ax.set_ylim([-0.1, 3.6])
(-0.1, 3.6)
%matplotlib inline
plt.figure(figsize=(15, 8), dpi=80, facecolor='w', edgecolor='k')
ax = plt.subplot(1, 1, 1)
ax.set_aspect('equal')
transform = matplotlib.transforms.Affine2D().rotate_deg(angle)
trans_data = transform + ax.transData
balkony = matplotlib.patches.Rectangle((0.0, 0.0), balcony_length, balcony_width, alpha=0.1, facecolor='black')
balkony.set_transform(trans_data)
ax.add_patch(balkony)
pv_module_1 = matplotlib.patches.Rectangle((0, 0), module_length, module_width, alpha=0.2, facecolor='yellow')
pv_module_1.set_transform(trans_data)
ax.add_patch(pv_module_1)
origin = np.array([0, 0])
yhat = np.array([0, 1])
x, y, arrow_length = 0.5, 0.5, 0.1
ax.annotate('N', xy=(x, y), xytext=(x, y-arrow_length), arrowprops=dict(facecolor='black', width=5, headwidth=15), ha='center', va='center', fontsize=20, xycoords=ax.transAxes)
ax.set_xlim([-1.2, 3.3])
ax.set_ylim([-0.1, 3.6])
(-0.1, 3.6)
lds = pd.Series(pd.date_range(start, freq=freq, periods=periods))
lds.head()
0 2023-08-11 12:30:00+02:00 1 2023-08-11 12:35:00+02:00 2 2023-08-11 12:40:00+02:00 3 2023-08-11 12:45:00+02:00 4 2023-08-11 12:50:00+02:00 dtype: datetime64[ns, Europe/Berlin]
def sun_elevation_and_azimuth(time, latitude, longitude, altitude, pressure=None, method='nrel_numpy', temperature=25.0, balkony_syste_rotate_right=0.0):
sun = pvlib.solarposition.get_solarposition(time, latitude, longitude, altitude=altitude, pressure=pressure, method=method, temperature=temperature) # , **kwargs
elevation = sun.elevation.values[0]
azimuth = sun.azimuth.values[0] + balkony_syste_rotate_right
return pd.Series([elevation, azimuth], index=['elevation', 'azimuth'])
def solar_normal_vector(elevation, azimuth):
r_elevation = scipy.spatial.transform.Rotation.from_euler('x', elevation, degrees=True)
r_azimuth = scipy.spatial.transform.Rotation.from_euler('z', -azimuth, degrees=True)
r_full = np.matrix(r_azimuth.as_matrix()) * np.matrix(r_elevation.as_matrix())
normal_vector = r_full * np.matrix([[0],[1],[0]])
normal_vector_array = np.squeeze(np.asarray(normal_vector))
return pd.Series(normal_vector_array, index=['x', 'y', 'z'])
def sun_shadow_point_projection(n, below=2.708 - 0.957):
scale_below = -below/n[2]
point_below = n * scale_below
return pd.Series(np.array(point_below), index=['px', 'py', 'pz'])
def solar_normal_vector_series(start_time, freq='5T', periods=10, balkony_syste_rotate_right=23.0, latitude=latitude, longitude=longitude, altitude=altitude, pressure=None, method='nrel_numpy', temperature=25):
start = pd.Timestamp(start_time)
if not(start.tzinfo is not None and start.tzinfo.utcoffset(start) is not None):
start = start.tz_localize(tzlocal.get_localzone())
lds = pd.Series(pd.date_range(start, freq=freq, periods=periods))
ldf = lds.apply(lambda t: sun_elevation_and_azimuth(t, latitude, longitude, altitude, pressure=pressure, method=method, temperature=temperature, balkony_syste_rotate_right=balkony_syste_rotate_right))
ldf.index = lds
ldf1 = ldf.apply(lambda row: solar_normal_vector(row['elevation'], row['azimuth']), axis=1)
ldf2 = ldf1.apply(lambda row: sun_shadow_point_projection(row[['x', 'y', 'z']]), axis=1)
ldf3 = (pd.Series(ldf.index).dt.hour.apply(lambda n: str(n).zfill(2)) + pd.Series(ldf.index).dt.minute.apply(lambda n: str(n).zfill(2))).to_frame('ts')
ldf3.index = lds
return pd.concat([ldf, ldf1, ldf2, ldf3], axis=1)
ldf = solar_normal_vector_series(start, freq=freq, periods=periods, balkony_syste_rotate_right=angle)
ldf.head()
elevation | azimuth | x | y | z | px | py | pz | ts | |
---|---|---|---|---|---|---|---|---|---|
2023-08-11 12:30:00+02:00 | 55.620088 | 181.756410 | -0.017308 | -0.564412 | 0.825312 | 0.036720 | 1.197470 | -1.751 | 1230 |
2023-08-11 12:35:00+02:00 | 55.907216 | 183.838478 | -0.037524 | -0.559277 | 0.828131 | 0.079342 | 1.182536 | -1.751 | 1235 |
2023-08-11 12:40:00+02:00 | 56.165516 | 185.947226 | -0.057691 | -0.553799 | 0.830650 | 0.121612 | 1.167402 | -1.751 | 1240 |
2023-08-11 12:45:00+02:00 | 56.394286 | 188.080364 | -0.077797 | -0.547980 | 0.832866 | 0.163560 | 1.152061 | -1.751 | 1245 |
2023-08-11 12:50:00+02:00 | 56.592883 | 190.235356 | -0.097834 | -0.541823 | 0.834779 | 0.205214 | 1.136505 | -1.751 | 1250 |
row = ldf[ldf['ts'] == '1300']
row
elevation | azimuth | x | y | z | px | py | pz | ts | |
---|---|---|---|---|---|---|---|---|---|
2023-08-11 13:00:00+02:00 | 56.897355 | 194.599614 | -0.137662 | -0.528506 | 0.837694 | 0.287749 | 1.104717 | -1.751 | 1300 |
%matplotlib inline
plt.figure(figsize=(15, 8), dpi=80, facecolor='w', edgecolor='k')
ax = plt.subplot(1, 1, 1)
ax.set_aspect('equal')
balkony = matplotlib.patches.Rectangle((0.0, 0.0), balcony_length, balcony_width, alpha=0.1, facecolor='black')
ax.add_patch(balkony)
pv_module_1 = matplotlib.patches.Rectangle((delta_x, delta_y), module_length, module_width, alpha=0.2, facecolor='yellow')
ax.add_patch(pv_module_1)
x1 = [delta_x+module_length/2.0, delta_x+module_length/2.0]
y1 = [0.0, module_width + 2*delta_y]
ax.plot(x1, y1, color='black')
for i,s in enumerate(['1250', '1300', '1345', '1420', '1540']):
row = ldf[ldf['ts'] == s]
rectangle_ = matplotlib.patches.Rectangle((row['px'], row['py']), balcony_length, balcony_width, alpha=0.2, facecolor='brown')
ax.add_patch(rectangle_)
ax.set_xlim([-1, 5])
ax.set_ylim([-0.5, 4])
(-0.5, 4.0)
location = pvlib.location.Location(latitude=latitude, longitude=longitude, tz='Europe/Berlin', altitude=350, name='Berghütte, 85863 Lansing')
location
Location: name: Berghütte, 85863 Lansing latitude: 48.1387900033257 longitude: 11.573658530854429 altitude: 350 tz: Europe/Berlin
up_to_date_cec_modules_df = None
def up_to_date_cec_modules():
global up_to_date_cec_modules_df
if up_to_date_cec_modules_df is not None:
return up_to_date_cec_modules_df
cec_modules_ = pvlib.pvsystem.retrieve_sam('CECMod')
cec_modules_filename = './config/CEC Modules.csv'
cec_modules = pd.read_csv(cec_modules_filename)
cec_modules = cec_modules.drop(axis=0, labels=[0, 1])
module_column_names = ['Technology', 'Bifacial', 'STC', 'PTC', 'A_c', 'Length', 'Width', 'N_s', 'I_sc_ref', 'V_oc_ref', 'I_mp_ref', 'V_mp_ref', 'alpha_sc', 'beta_oc', 'T_NOCT', 'a_ref', 'I_L_ref', 'I_o_ref',
'R_s', 'R_sh_ref', 'Adjust', 'gamma_r', 'BIPV', 'Version', 'Date']
module_column_types = [str, int, float, float, float, float, float, int, float, float, float, float, float, float, float, float, float, float, float, float, float, float, str, str, str]
for i, c in enumerate(module_column_names):
cec_modules[c] = cec_modules[c].astype(module_column_types[i])
# cec_module_names = list(cec_modules['Name'])
cec_modules = cec_modules.set_index('Name')
cec_modules = cec_modules[cec_modules_.index]
cec_modules = cec_modules.T
up_to_date_cec_modules_df = cec_modules
return cec_modules
Where and how to find 'CEC Inverters.csv' and 'CEC Modules.csv':
up_to_date_cec_inverters_df = None
def up_to_date_cec_inverters():
global up_to_date_cec_inverters_df
if up_to_date_cec_inverters_df is not None:
return up_to_date_cec_inverters_df
cec_inverters_ = pvlib.pvsystem.retrieve_sam('CECInverter')
cec_inverters_filename = './config/CEC Inverters.csv'
cec_inverters = pd.read_csv(cec_inverters_filename) # skiprows=2, index_col=0,
cec_inverters = cec_inverters.drop(axis=0, labels=[0, 1])
inverter_column_names = ['Name', 'Vac', 'Pso', 'Paco', 'Pdco', 'Vdco', 'C0', 'C1', 'C2', 'C3', 'Pnt', 'Vdcmax', 'Idcmax', 'Mppt_low', 'Mppt_high', 'CEC_Date', 'CEC_hybrid']
inverter_column_types = [str, float, float, float, float, float, float, float, float, float, float, float, float, float, float, str, str]
for i, c in enumerate(inverter_column_names):
cec_inverters[c] = cec_inverters[c].astype(inverter_column_types[i])
cec_inverters['CEC_Type'] = 'Utility Interactive' # array(['Utility Interactive', 'Grid Support'], dtype=object)
cec_inverters = cec_inverters.set_index('Name')
cec_inverters = cec_inverters[cec_inverters_.index]
cec_inverters = cec_inverters.T
up_to_date_cec_inverters_df = cec_inverters
return cec_inverters
def cec_system(surface_tilt_south = 0.0):
south_modules = 1
cec_modules = up_to_date_cec_modules()
cec_inverters = up_to_date_cec_inverters()
module = cec_modules['Jinko Solar Co. Ltd JKM375M-72']
inverter = cec_inverters['Hoymiles Power Electronics Inc : HMS-600-2D-NA [240V]']
temperature_parameters = pvlib.temperature.TEMPERATURE_MODEL_PARAMETERS['sapm']['open_rack_glass_glass']
# Süd: 157° West: 247°
mount_south = pvlib.pvsystem.Array(pvlib.pvsystem.FixedMount(surface_tilt=surface_tilt_south, surface_azimuth=157), module_parameters=module, modules_per_string=south_modules, strings=1, temperature_model_parameters=temperature_parameters)
system = pvlib.pvsystem.PVSystem(arrays=[mount_south], inverter_parameters=inverter)
model_chain = pvlib.modelchain.ModelChain(system, location, aoi_model="no_loss") # aoi_model='physical', spectral_model='first_solar'
return model_chain
model_chain = cec_system()
pd.Timestamp(start.date()) + pd.Timedelta('6H')
Timestamp('2023-08-11 06:00:00')
times = pd.date_range(start=(pd.Timestamp(start.date()) + pd.Timedelta('6H')), end=(pd.Timestamp(start.date()) + pd.Timedelta('20H')), freq='T', tz=location.tz)
pvgis_minute = location.get_clearsky(times)
pvgis_minute['precipitable_water'] = 1.42
pvgis_minute.head()
ghi | dni | dhi | precipitable_water | |
---|---|---|---|---|
2023-08-11 06:00:00+02:00 | 0.0 | 0.0 | 0.0 | 1.42 |
2023-08-11 06:01:00+02:00 | 0.0 | 0.0 | 0.0 | 1.42 |
2023-08-11 06:02:00+02:00 | 0.0 | 0.0 | 0.0 | 1.42 |
2023-08-11 06:03:00+02:00 | 0.0 | 0.0 | 0.0 | 1.42 |
2023-08-11 06:04:00+02:00 | 0.0 | 0.0 | 0.0 | 1.42 |
plt.figure(figsize=(32, 8), dpi=80, facecolor='w', edgecolor='k')
major_ticks_y = np.arange(0, 301, 25)
ax = plt.subplot(1, 1, 1)
ax.set_yticks(major_ticks_y)
minlocator = matplotlib.dates.MinuteLocator(byminute=list(range(0,60,20)), tz=location.tz)
minlocator.MAXTICKS = 40000
majorFmt = matplotlib.dates.DateFormatter('%H:%M:%S', tz=location.tz)
ax.xaxis.set_major_locator(minlocator)
ax.xaxis.set_major_formatter(majorFmt)
plt.setp(ax.xaxis.get_majorticklabels(), rotation=90)
ax.grid(True)
ldf = model_chain.run_model(pvgis_minute).results.dc.p_mp.to_frame('dc')
ax.plot(ldf.index, ldf['dc'])
ax.set_xlim([ldf.index[0],ldf.index[-1]])
plt.show()
ldf['dc'].max()
270.5100350769249
dc_production = (ldf * 1/60.0).sum()
dc_production
dc 2274.546351 dtype: float64
# compared to 1566 Wh produced
real_dc_production = 1566.0
loss_absolute = dc_production - real_dc_production
loss_absolute
dc 708.546351 dtype: float64
loss_percent = loss_absolute / dc_production
loss_percent
dc 0.311511 dtype: float64
increase_potential = dc_production / real_dc_production
round(increase_potential,2)
dc 1.45 dtype: float64