Schattingen van jaaropbrengsten van zonnepanelen zijn gemakkelijk op internet te vinden, maar hoe ziet de opwek per maand eruit voor een specifieke opstelling, per dag of zelfs per uur? Python heeft een pakket om zonnepanelen te simuleren: pvlib. In deze post lees je hoe het werkt!

In dit voorbeeld worden zonnepanelen gesimuleerd die op het zuiden liggen onder een hellingshoek van 30°, en een vermogen van 1000 wattpiek (ofwel 1 kWp). Het simuleren met pvlib bestaat uit vier stappen:

  • Weerdata ophalen van het KNMI
  • Diffuse en directe straling schatten aan de hand van de globale straling
  • De zonnepanelen-installatie definiëren
  • Opwek simuleren

Als pvlib eenmaal is geïnstalleerd, moeten eerst de benodigde pakketten worden ingeladen: 

import pandas as pd
import matplotlib.pyplot as plt
import datetime
import requests

from pvlib import pvsystem, modelchain, location, irradiance
from pvlib.solarposition import get_solarposition

Weerdata ophalen van het KNMI

Om de opwek van zonnepanelen te simuleren, heb je gegevens over de zon nodig: hoeveel zonnestraling is er op een gegeven moment? Door KNMI-weerdata hiervoor te gebruiken, kun je de opwek in een historische periode simuleren.

Met onderstaande functie worden de juiste variabelen – zonnestraling en temperatuur – opgehaald bij het KNMI. Vervolgens worden deze variabelen in de eenheden gezet zoals pvlib ze verwacht. Als tijdstip bij elk datapunt wordt de tijd halverwege het gemeten uur genomen. 

def get_hourly_weather_data_for_pvlib(stations, start_date, end_date, timezone = 'UTC'):
    '''
    Function to get hourly weather variables T (temperature) and
    Q (global radiation) from KNMI 
    
    Args: 
        stations   (str): NMI-stations separated by ':' 
        start_date (str): start date, format yyyymmdd
        end_date   (str): end date (included), format yyyymmdd
        timezone   (str, optional): timezone

    Returns:
        df: DataFrame with DateTime-index, columns T (temp), Q (global radiation) 
    '''

    url = 'https://www.daggegevens.knmi.nl/klimatologie/uurgegevens'

    data = {
        'start': start_date,
        'end': end_date,
        'vars': 'Q:P:T',
        'stns': stations,
        'fmt': 'json'
        }

    response = requests.post(url, data = data)    
    weather_df = pd.DataFrame(response.json())

    # correct units
    weather_df['T'] = weather_df['T'] / 10          # is in 0.1 degrees C, to degrees C
    weather_df['Q'] = weather_df['Q'] * (1 / 0.36)  # is in J/m2, to W / m2
    weather_df['P'] = weather_df['P'] * 10          # is in 0.1 hPa, to Pa
    
    # create date_time index, convert timezone
    weather_df['hour'] = weather_df['hour'] - 1     # is from 1-24, to 0-23
    weather_df['date_time'] = pd.to_datetime(weather_df['date']) +
                              pd.to_timedelta(weather_df['hour'].astype(int), unit='h')
    weather_df.index = weather_df.date_time
    weather_df = weather_df.drop(['date', 'hour', 'date_time'], axis = 1)
    weather_df.index = weather_df.index.tz_convert(timezone)

    # shift date_time by 30 minutes, 'average time' during that hour
    weather_df.index = weather_df.index.shift(freq="30min")

    return weather_df

In het voorbeeld wordt data van het weerstation in De Bilt gebruikt om de opwek in 2023 te simuleren. De lengte- en breedtegraad van De Bilt worden gebruikt om straks de positie van de zon te kunnen bepalen. 

start_date = '20230101'
end_date = '20231231'

timezone = 'Europe/Amsterdam'

# De Bilt 
station = '260'
lat = 52.1092717
lon = 5.1809676

# Function get_hourly_weather_data_for_pvlib defined in full code overview below
weather_df = get_hourly_weather_data_for_pvlib(station, start_date, end_date, timezone)

Van globale straling naar directe en diffuse straling

Nu weten we de globale straling per uur gedurende het jaar, maar we hebben nog meer nodig om zonnepanelen te simuleren: er moet een schatting gemaakt worden van de zogeheten DNI (direct normal irradiance) en DHI (diffuse horizontal irradiance). Er zijn in pvlib verschillende methodes om aan de hand van de globale straling (GHI, ofwel global horizontal irradiance) een schatting te maken van DNI en DHI. Hieronder is gekozen voor de methode genaamd Erbs. Deze functie maakt gebruik van de stand van de zon in De Bilt, die wordt opgehaald met de get_solarposition functie van pvlib. 

# Get solar position for date_times
solpos_df = get_solarposition(
    weather_df.index, latitude = lat,
    longitude = lon, altitude = 0,
    temperature = weather_df['T'])
solpos_df.index = weather_df.index

# Use method 'Erbs' to go from GHI to DNI and DHI. There are other methods available.
irradiance_df = irradiance.erbs(weather_df['Q'], solpos_df['zenith'], weather_df.index)
irradiance_df['ghi'] = weather_df['Q']

Zonnepanelen definiëren

De dataframe irradiance_df bevat nu de juiste kolommen om opwek te simuleren, dus wordt het tijd om de zonnepanelen zelf te definiëren. Hoewel je met pvlib allerlei parameters van de zonnepanelen aan kan passen, zijn er vier vooral van belang: de oriëntatie van de panelen, de hellingshoek, het totaal aantal Wattpiek (Wp) van de panelen en het vermogen van de omvormer.

De panelen in dit voorbeeld liggen pal op het zuiden onder een hellingshoek van 30°. In totaal ligt er 1000 Wp aan panelen, en de omvormer in het voorbeeld is ook 1000 W. 

# Define solar panels
orientation = 180  # Azimuth, 90 is East, 180 is South, 270 is West
tilt = 30         # 0 is flat on the ground, 90 is vertical
Wp_panels = 1000  # total Wp of the solar panels
W_inverter = 1000 # power (W) of the inverter

array = pvsystem.Array(pvsystem.FixedMount(tilt, orientation),
                       module_parameters=dict(pdc0 = Wp_panels, gamma_pdc = -0.004),
                       temperature_model_parameters = {'a': -3.56, 'b': -.0750, 'deltaT': 3})


loc = location.Location(lat, lon, tz = timezone)
system = pvsystem.PVSystem(arrays = [array], inverter_parameters = dict(pdc0 = W_inverter))
mc = modelchain.ModelChain(system, loc, aoi_model = 'physical',
                           spectral_model = 'no_loss')

Eerst wordt de array gedefinieerd: de zonnepanelen zelf. De zonnepanelen liggen op een vaste plek, vandaar dat er voor FixedMount() gekozen is. Er moeten daarnaast nog een aantal parameters worden meegegeven, waaronder temperature parameters. Zonnepanelen worden minder efficiënt als ze te warm worden, en hoe warm ze worden hangt onder andere af van het soort panelen (glas-glas of glas-folie) en of de zonnepanelen in een open rack zitten, of in het dak. Zie de documentatie voor de opties, hier zijn de waardes gekozen voor glas-glas open rack. De parameter gamma_pdc geeft aan hoe de vermogensoutput van een zonnepaneel verandert met de temperatuur, en zit volgens pvlib over het algemeen tussen de -0.002 and -0.005 per graad Celsius.

Het totale systeem van de zonnepanelen bestaat uit de array en de omvormer (inverter). Dit systeem kan worden meegegeven in de ModelChain, waar ook de locatie wordt meegegeven.

Opwek simuleren

Nu kan het simuleren van de opwek dan echt beginnen. Dit kan door stralingsdata aan het model te voeren. Eerst proberen we het voor één dag. 

date = datetime.date(2023, 6, 10)
mc.run_model(irradiance_df[irradiance_df.index.date == date])
mc.results.ac.plot()
plt.title(f'pvlib: opwek gesimuleerde zonnepanelen (1 kWp) op {date}', size  = 12)
plt.ylabel('Opwek (W)', size = 12)
plt.xlabel('Tijdstip', size = 12)
plt.show()

# Compute total production in kWh
pvlib_df = mc.results.ac.reset_index()
pvlib_df['kWh'] = pvlib_df['p_mp'] / 1000  
production_kWh = round(pvlib_df['kWh'].sum(),2)
print(f'On {date} the solar panels produced {production_kWh} kWh')

De gekozen datum (10 juni 2023) is een dag die bijna perfect zonnig was. De zonnepanelen zouden 7,21 kWh hebben opgewekt op deze datum. De figuur die je met de plot-functie van pvlib kan maken, ziet er zo uit:

Opwek zonnepanelen met pvlib

Op dezelfde manier kan pvlib het hele jaar simuleren. Hieruit blijkt dat deze zonnepanelen in 2023 1107,19 kWh zouden hebben opgewekt:

mc.run_model(irradiance_df)

# Compute total production in kWh
pvlib_df = mc.results.ac.reset_index()
pvlib_df['kWh'] = pvlib_df['p_mp'] / 1000  
production_kWh = round(pvlib_df['kWh'].sum(),2)
print(f'In 2023, the solar panels produced {production_kWh} kWh')

Door verschillende opstellingen te proberen, kun je zo dus kijken wat de verschillen zijn in jaarverbruik als de panelen niet op het zuiden liggen, of onder een andere hoek.

Zonnepanelen simuleren zonder te programmeren

Zonnepanelen simuleren maar niet programmeren? Dat kan met PVGIS, een gratis webapplicatie gemaakt door het Joint Research Center van de EU. Maandelijkse, dagelijkse en uurlijkse opwek is te downloaden als csv. Wel gaat de weerdata hier maar tot en met 2020. Deze tool is overigens ook via pvlib aan te roepen.

Typisch Nederlands weer?

Een andere optie dan historische weerdata van het KNMI gebruiken, is om gegevens van een TMY te gebruiken, ofwel een typical meteorological year: een gemiddeld weerjaar voor een gegeven locatie. Op de site van PVGIS is een bestand met gegevens van een TMY van een zelfgekozen locatie te downloaden. Hiermee heb je er geen last van dat de door jou gekozen periode misschien extreem zonnig was, of juist bovengemiddeld bewolkt.

Volledige code zonnepanelen simuleren met Python’s pvlib

import pandas as pd
import matplotlib.pyplot as plt
import datetime
import requests

from pvlib import pvsystem, modelchain, location, irradiance
from pvlib.solarposition import get_solarposition


def get_hourly_weather_data_for_pvlib(stations, start_date, end_date, timezone = 'UTC'):
    '''
    Function to get hourly weather variables T (temperature) and
    Q (global radiation) from KNMI 
    
    Args: 
        stations   (str): NMI-stations separated by ':' 
        start_date (str): start date, format yyyymmdd
        end_date   (str): end date (included), format yyyymmdd
        timezone   (str, optional): timezone

    Returns:
        df: DataFrame with DateTime-index, columns T (temp), Q (global radiation) 
    '''

    url = 'https://www.daggegevens.knmi.nl/klimatologie/uurgegevens'

    data = {
        'start': start_date,
        'end': end_date,
        'vars': 'Q:P:T',
        'stns': stations,
        'fmt': 'json'
        }

    response = requests.post(url, data = data)    
    weather_df = pd.DataFrame(response.json())

    # correct units
    weather_df['T'] = weather_df['T'] / 10          # is in 0.1 degrees C, to degrees C
    weather_df['Q'] = weather_df['Q'] * (1 / 0.36)  # is in J/m2, to W / m2
    weather_df['P'] = weather_df['P'] * 10          # is in 0.1 hPa, to Pa
    
    # create date_time index, convert timezone
    weather_df['hour'] = weather_df['hour'] - 1     # is from 1-24, to 0-23
    weather_df['date_time'] = pd.to_datetime(weather_df['date']) +
                              pd.to_timedelta(weather_df['hour'].astype(int), unit='h')
    weather_df.index = weather_df.date_time
    weather_df = weather_df.drop(['date', 'hour', 'date_time'], axis = 1)
    weather_df.index = weather_df.index.tz_convert(timezone)

    # shift date_time by 30 minutes, 'average time' during that hour
    weather_df.index = weather_df.index.shift(freq="30min")

    return weather_df


start_date = '20230101'
end_date = '20231231'

timezone = 'Europe/Amsterdam'

# De Bilt 
station = '260'
lat = 52.1092717
lon = 5.1809676

weather_df = get_hourly_weather_data_for_pvlib(station, start_date, end_date, timezone)

# Get solar position for date_times
solpos_df = get_solarposition(
    weather_df.index, latitude = lat,
    longitude = lon, altitude = 0,
    temperature = weather_df['T'])
solpos_df.index = weather_df.index

# Use method 'Erbs' to go from GHI to DNI and DHI. There are other methods available.
irradiance_df = irradiance.erbs(weather_df['Q'], solpos_df['zenith'], weather_df.index)
irradiance_df['ghi'] = weather_df['Q']


# Define solar panels
orientation = 180  # Azimuth, 90 is East, 180 is South, 270 is West
tilt = 30         # 0 is flat on the ground, 90 is vertical
Wp_panels = 1000  # total Wp of the solar panels
W_inverter = 1000 # power (W) of the inverter


array = pvsystem.Array(pvsystem.FixedMount(tilt, orientation),
                       module_parameters = dict(pdc0 = Wp_panels, gamma_pdc = -0.004),
                       temperature_model_parameters = {'a': -3.56, 'b': -.0750, 'deltaT': 3})

loc = location.Location(lat, lon, tz = timezone)
system = pvsystem.PVSystem(arrays = [array], inverter_parameters = dict(pdc0 = W_inverter))
mc = modelchain.ModelChain(system, loc, aoi_model = 'physical',
                           spectral_model = 'no_loss')


# Example 1 day
date = datetime.date(2023, 6, 10)
mc.run_model(irradiance_df[irradiance_df.index.date == date])
mc.results.ac.plot()
plt.title(f'pvlib: opwek gesimuleerde zonnepanelen (1 kWp) op {date}', size  = 12)
plt.ylabel('Opwek (W)', size = 12)
plt.xlabel('Tijdstip', size = 12)
plt.show()

# Compute total production in kWh
pvlib_df = mc.results.ac.reset_index()
pvlib_df['kWh'] = pvlib_df['p_mp'] / 1000  
production_kWh = round(pvlib_df['kWh'].sum(),2)
print(f'On {date} the solar panels produced {production_kWh} kWh')


# Example whole year
mc.run_model(irradiance_df)

# Compute total production in kWh
pvlib_df = mc.results.ac.reset_index()
pvlib_df['kWh'] = pvlib_df['p_mp'] / 1000  
production_kWh = round(pvlib_df['kWh'].sum(),2)
print(f'In 2023, the solar panels produced {production_kWh} kWh')