Ziel

Im Tutorial 2 von https://pv4ing.ch sind drei Übungen vorgesehen. Die erste Übung erstellt eine Zeitreihe mit fiktiven Daten. Die zweite Übung verwendet Messdaten als Zeitreihe und die dritte Übung führt eine Berechnung der Jahresstrahlungsenergie durch, bei unterschiedlichen Ausrichtungen.

Übung 1: Einstieg

Zuerst machen wir eine einfache Zeitreihenrechnung.

Schritt 1 - Synthetische Messwerte für 1 Tag

Im folgenden Beispiel, nehmen wir keine realen Messwerte, sondern "bauen" uns selbst die Messwerte über die Sinusfunktion. Zuerst erstellen wir einen Array (Zahlenreihe) für die Zeit. Wir nehmen an, dass die Sonne um 6:00 aufgeht und um 18:00 untergeht. Über Mittag erreicht die Strahlung einen Wert von 650 W/m2.

In [1]:
# Bibliotheken Import
import numpy as np
import matplotlib.pyplot as plt
import datetime as dt

# Funktionsdefinition
cos  = lambda arg : np.cos(np.deg2rad(arg))
sin  = lambda arg : np.sin(np.deg2rad(arg))
acos = lambda arg : np.rad2deg(np.arccos(arg))
asin = lambda arg : np.rad2deg(np.arcsin(arg))

t = np.linspace(6,18,1000)
H = 650 * sin((t-6)/(18-6)*180)

Wir haben ein Strahlungsarray H generiert, eine Zahlenreihe mit Strahlungswerte in W/m2. Beachte das die Funktion sin() oben mit dem lambda-Operator definiert wurde und in Grad angegeben wird.

Überprüfe ob der Zeitarray aus lauter Spaltenwerte (Zeile mit vielen Zahlen) oder aus lauter Zeilenwerte (Spalte mit vielen Zahlen) besteht. Spalte = column. Zeile = row.

In [2]:
t.shape
Out[2]:
(1000,)

Stelle den Strahlungsverlauf als Liniendiagramm mit der Linienfarbe grün (g=green) dar.

In [3]:
plt.figure(2, figsize=(8,4))
plt.plot(t,H, 'g')
plt.xlabel('Zeit [h]')
plt.ylabel('Strahlung [W/m2]')
Out[3]:
Text(0, 0.5, 'Strahlung [W/m2]')

Nun können wir die Strahlungsenergie $W_H$ für diesen einen Tag berechnen. Beachte das die Berechung der "Energie = Leistung mal Zeit" ist. Hier ist der Zeitschritt "deltaT" relevant.

Die Jahresstrahlungsenergie geben wir als String aus in der Form: "Strahlungsenergie 5 kWh auf 1.0 m2 Fläche in einem Tag"

Es wird ein Sting erstellt indem das Rechenergebnis wH und die Fläche A eingebaut wird. Die Zahlen werden mit {} im Text angegeben. Dieser {}-Bereich wird ersetzt mit Zahlenwerten welche im .format() angegeben werden. In unserem Fall sind es die beiden Werte A und wH. In der ersten Klammer {1:3.0f} geben wir an, dass wir das zweite Element mit Index=1 wollen. Das erste Element hat der Index=0. Dann folgt ein ":" und die Formatierung:

3.0f heisst eine Float-Zahl mit drei Stellen und 0 Nachkommastellen.

In [4]:
A = 1 # Fläche 1 m2
deltaT = t[1] - t[0]
wH = np.sum(H)*deltaT*A/1000
print('Strahlungsenergie {1:3.0f} kWh auf {0:1.1f} m2 Fläche in einem Tag'.format(A,wH))
Strahlungsenergie   5 kWh auf 1.0 m2 Fläche in einem Tag

Bildhafte Darstellung der Funktionsweise von .format. textformat

Schritt 2 - Synthetische Messwerte für 7 Tag

Nun wollen wir die Strahlungsenergie von 7 Tagen berechnen, von Stunde 0 bis Stunde 7 x 24. Die Sinusfunktion liefert uns für die Nacht (18:00 bis 6:00) negative Werte. Wir korrigieren die negativen Strahlungswerte zu 0 indem wir die Indexierung von Python verwenden.

In [5]:
H[H<0] = 0

Der Ausdruck "H<0" generiert ein Selektor-Array der Länge von H, bei dem die einzelnen Positionen mit "1" oder "True" gekennzeichnet sind, bei welchen "Wert < 0" zutrifft. Der Ausdruck H[Selektor-Array] = 0 weisst allen selektierten Positionen eine 0 zu.

Denkbar ist auch "M[N<1] = 2". M und N können unterschiedliche Inhalte haben, müssen aber gleich lange Arrays sein. Im Array M werden alle Stellen auf 2 gesetzt, bei denen der Wert im Array "N < 1" zutriff.

Nun haben wir folgende Punkte gelernt:

  1. Erstellen eines Arrays.
  2. Liniendiagramm inklusiv Beschriftung.
  3. Berechnen eines Flächenintegrals aus diskreten Werten, d.h. Aufsummieren von Werten.
  4. Ergebnisse als formatierbare Zahlen zusammen mit Text ausgeben.
  5. Erstellen eines Array über mehrere Tag mit der Verwendung der Indexierung

Übung 2: Berechnung Strahlungsenergie von mehreren Tagen mit Messdaten

Nun verwenden wir nicht "künstliche" Strahlungswerte, generiert durch die Sinusfunktion, sondern reale Messwert. Hierzu benötigen wir eine Datendatei. Wir erstellen einen Plot über 7 Tage mit der Globalstrahlung "hGlob", der Diffusstrahlung "hDif" und der Direktstrahlung "hDir".

Bei Energiesimulationen ist es empfehlenswert als Kommentar die Einheit anzugeben, sodass verständlich ist um was es sich handelt, z.B. [grad] oder [radiant] oder [W] oder [kWh].

Verwende für diese Übung ein neues py-File und importiere wieder die benötigten Bibliotheken (numpy, matplotlib, datetime).

Datenimport mit Pandas

Zum Importieren verwenden wir die Bibliothek Pandas. Diese zeichnet sich durch gute Verarbeitung von Daten aus, vorallem Zeitreihen. Das Einlesen geht mittels read_csv, dabei erhalten wir ein DataFrame (df). Dies ist vergleichbar mit einem Excel-Tabellenblatt mit einem Header, d.h. Überschriftzeile und Links einem laufenden Index.

Als Übung verwenden wir Daten von PVGIS über das Tool "TMY". Wir wählen einen Ort und laden die Daten als csv-Datei. pvgis

Sieh dir die csv-Datei mit dem Texteditor an. Sie enthält einige Information zu Beginn und am Ende. Bei Einlesen sollen dies Informationen ignoriert werden. Dies wird mit skiprows und skipfooter angegeben. Diese Funktionen sind mit engine = 'python' verfügbar. Standardmässig verwendet read_csv die 'c'-engine, da diese schneller ist, jedoch nicht so mächtig.

In [6]:
import pandas as pd

df = pd.read_csv('pvgis.csv', delimiter=',', decimal='.', skiprows = 16, skipfooter = 11, header = 0, engine = 'python')
df['tutc'] = pd.to_datetime(df['time(UTC)'], format='%Y%m%d:%H%M' )

Pandas erkennt automatisch, welche Daten "strings" sind und welche Daten "Zahlen", d.h. float sind. Pandas erkennt aber nicht das Datumsformat. Hier helfen wir, indem wir das Format angeben und die Zeitwerte als maschinenlesbares Format in eine neue Spalte tutc speichern.

Nach dem Import, schau dir die importen Daten als DataFrame (df) an, mit

  • print(df) oder
  • print(df.shape) oder
  • list(df) oder
  • describe(df).

Zeitreihenrechnung

Zum Rechnen verwenden wir Numpy-Array, d.h. wir nehmen die Daten aus dem DataFrame (df) in dem wir über den Spaltennamen, in eckiger Klammer, auf die Spalte zugreifen und die Werte (values) auslesen:

In [7]:
hGlo  = df['G(h)'].values    # [W/m2] 
hDif  = df['Gd(h)'].values   # [W/m2]
hDNI  = df['Gb(n)'].values   # [W/m2]
tAmb  = df['T2m'].values     # [°C]

Für die Sonnenstandsberechnung benötigen wir einen Array mit laufenden Stunden "lfStd", beginnend mit Stunde 0 bei Jahresanfang. Die PVGIS-Daten über TMY (Typical Meteorological Year) wurden aus Monatsmessdaten von unterschiedlichen Jahren zusammengestellt. Wir verwenden die Daten für eine Jahressimulation, weshalb wir einen neuen Zeit-Array erstellen, mit der selben Länge wie die Daten.

In [8]:
lfStd = np.linspace(0,8760, 8760)
deltaT = 1 # Zeitintervall 1 Stunde

Schritt 1

Es soll die Globalstrahlung in grün, Direktstrahlung in rot und Diffusstrahlung in blau dargestellt werden. Weiter soll in der Plot-Figure eine Legende mit den Kurvennamen dargestellt werden. Als Ausschnitt wollen wir die ersten 7 Tage betrachten, d.h. in der Plotfunktion grenzen wir den Array ein durch lfStd[a:e], wobei a für den Anfangsindex und e für den Endindex steht.

Schritt 2

Wenn wir den ersten Tag im Plot ansehen, so fällt auf, dass negative Strahlungswerte vorkommen. Dies entspricht den Messergebnissen aber nicht der Realität. Deshalb sollen die Strahlungsdaten bereinigt werden.  

In [9]:
hDir = hGlo - hDif

hDir[hDir<0] = 0
hDif[hDif<0] = 0
hGlo[hGlo<0] = 0

a=int(0)           # Anfangstag
e=int(7*24)        # Endtag
t = lfStd[a:e]

plt.figure(3, figsize=(8,4))# Grösse des Plots (figsize) in Zoll
plt.plot(t,hGlo[a:e],'g', label ='Globalstrahlung')
plt.plot(t,hDir[a:e],'r', label ='Direktstrahlung')
plt.plot(t,hDif[a:e],'b', label ='Diffusstrahlung')
plt.xlabel('Zeit [h]');
plt.ylabel('Strahlung [W/m2]');
plt.legend(loc="upper left")
Out[9]:
<matplotlib.legend.Legend at 0x20456cec970>

Schritt 3

Nun berechnen wir die Strahlungsenergie über das gesamte Jahr. Die Ausgabe der Zahlen soll ohne Nachkommastellen erfolgen, deshalb wird die Zahl als formatierter String in den Gesamtstring eingebettet mit der bekannten Methode ".format".

In [10]:
wGlo = np.sum(hGlo)*deltaT/1000 # [kWh/m2] Jahresstrahlungsenergie Globalstrahlung
wDir = np.sum(hDir)*deltaT/1000 # [kWh/m2] Jahresstrahlungsenergie Direktstrahlung
wDif = np.sum(hDif)*deltaT/1000 # [kWh/m2] Jahresstrahlungsenergie Diffusstrahlung
print('Global {0:1.0f} kWh/m2 entspricht {1:3.0f} %'.format(wGlo, wGlo/wGlo*100))
print('Direkt {0:1.0f} kWh/m2 entspricht {1:3.0f} %'.format(wDir, wDir/wGlo*100))
print('Diffus {0:1.0f} kWh/m2 entspricht {1:3.0f} %'.format(wDif, wDif/wGlo*100))
Global 1237 kWh/m2 entspricht 100 %
Direkt 692 kWh/m2 entspricht  56 %
Diffus 545 kWh/m2 entspricht  44 %

Zur Kontrolle: vergleiche die Jahresglobalstrahlungsenergie für diesen Standort mit den Werten aus PVGIS.

Nun haben wir folgende Punkte gelernt:

  1. Umweltdaten in einem maschinenlesbaren Format zu organisieren über PVGIS und TMY
  2. Importieren von Messdaten mit Pandas.
  3. Bereinigen eines Array (z.B. Nullsetzten bei negativen Werten) mit der Indexierung.

Übung 3: Strahlung auf eine geneigte Fläche

Bisher betrachteten wir die Strahlung auf die Horizontale. Horizontal weil die Messgeräte horizontal aufgestellt sind. Über die Aufteilung der Globalstrahlung in Diffus- und Direktstrahlung kann die Strahlung auf eine beliebige Fläche berechnet werden.

Schritt 1

Berechne für einen beliebigen Ort die Jahresstrahlungsenergie bei unterschiedlichen Neigungen. Verwende PVGIS als Datenquelle für Strahlungsdaten. Berechne mit dem 3K-Modell (Dreikomponentenmodell) die Strahlung auf eine beliebig ausgerichtete Fläche.Das Dreikomponentenmodell und die Sonnenstandsberechnung ist im Manuskript "Photvoltaik für Ingenieure" beschrieben.

Ziel ist ein Liniendiagramm mit der Jahresstrahlungsenergie bei einer Modulausrichtung nach Süden (Azimut 0°). Auf der horizontalen x-Achse ist die Modulneigung und auf der vertikalen y-Achse die Jahresstrahlungsenergie.

Aufbau des Code-Beispiels:

Mit dem 3K-Modell wird die Strahlungsleistung übers gesamte Jahr berechnet. Am Ende wird die Jahresenergie (wH) mit deltaT multipliziert und aufsummieren. Wir erhalten die Jahresstrahlungsenergie

Um diese Code-Zeilen wird eine for-Schleife gelegt um die Neigung zu variieren.

Weitere Hinweise zum Code:

  • Bei der Berechnung des Direktstrahlungsanteils ist eine Division von sin(h) notwendig. Wenn h = 0 ist, führt dies zu Problemen in der Berechnung, weshalb die Werte auf eine Sonnenhöhe von $5^\circ$ limitiert werden. Implementierung sieht wie folgt aus, anstatt sin(h) wird np.maximum(sin(h),sin(5)) verwendet.
  • Bei der Python-Implementierung ist darauf zu achten das "np.maximum()" verwendet wird, anstatt "max"! Diese Funktionen sind nicht identisch. Überprüfe den Unterschied. np.maximum ermittelt elementweise das Maximum bei zwei Array's.
  • Bei der Berechnung der Direktstrahlungskomponente wird der Zwischenwinkel $\Theta$ benötigt. Ist dieser über 90° heisst dies das die Direktstrahlung hinter dem Modul ist. Deshalb müssen diese Werte korrigiert werden.

    theta[theta>90] = 90

    Bei genauerer Betrachtung der Berechnung ist nicht $\Theta$ in der Berechnung notwendig, sondern $cos(\Theta)$. Die Korrektur sieht dann wie folgt aus:

    cosTheta[cosTheta<0] = 0

Schritt 2

Zur Kontrolle:

  • Stelle den Sonnenstand dar (Azimut und Sonnenhöhe)
  • Vergleiche die Ergebnisse mit PVGIS. Beachte das die Jahresstrahlungsenergie in kWh/m2 angegeben wird. In PVGIS ist dies "H" nicht "E". Trage die PVGIS-Werte ebenfalls ins Diagramm ein.

Schritt 3

Aufgabe Optimale Ausrichtung: Variiere nun ein weiterer Parameter. Neben der Neigung soll die Ausrichtung, d.h. der Azimut variiert werden. Auch hier wieder Jahresberechnung der Jahresstrahlungsenergie. Das Ergebnis ist ein 3D-Plot oder eine Höhenliniendarstellung auch Contour-Plot genannt.

Zusammenfassung Tutorial 2

Übung 1

Energie ist Leistung mal Zeit. Leistungswerte werden summiert und mit DeltaT multipliziert.

Übung 2

Importieren von Daten und Darstellen.

Übung 3

Die Berechnung der Jahresenergie erfolgt vektoriell.

Zwei ineinander verschachtelte for-Schleifen! Es werden zwei Parameter variiert, deshalb x-, y-, z-Achse. Für jeden Datenpunkt erfolgt eine Jahressimulationen. Im Beispiel wurde die Neigungen und die Ausrichtungen variiert deshalb zwei for-Schleifen.

Weitere Informationen zu Energiesystemberechnung

Bisher verwendeten wir die Umweltdaten Strahlung und Tempertaur. Für Berechnungen zu Eigenverbrauchsquote, Autarkie, Warmwassererwärmung benötigen wir Lastprofiele. Das Lastprofil ist syntetisch erstellt und sieht wie folgt aus:

Wir verwenden die csv-Datei auf github https://github.com/markstaler/pv4ing. Mit "save_as" kann die Datei gespeicher werden. Achtung, wenn eine csv-Datei angeklickt wird, so wird diese standardmässig mit Excel geöffnet, dabei formatiert Excel die Zahlen, Datum, Uhrzeit um, mit den Einstellungen des Betriebssystems. Dies geht in 9 von 10 Fällen gut, aber nicht immer. Deshalb die Daten nicht speichern, beim Öffnen mit Excel. Besser die Datei mit dem Editor/Notepad öffnen.

Im DataFrame ist das Datum und die Uhrzeit angegeben, jedoch als Text-String und für pandas nicht verständlich. Dieser Text-String wandeln wir in ein maschinenlesbares Format um, hierbei geben wir das Format des Datums und der Uhrzeit an.

Anschliessend legen wir diese Spalte mit Datum und Uhrzeit als Index fest (set_index) und die ursprüngliche Spalte "Time" wird gelöscht. Mit der Festlegung von datetime als Index weiss pandas dass es sich um eine Zeitreihe handelt und kann diese, sofern gewünscht in eine beliebige Auflösung umrechnen mit resample.

Im untenstehenden Beispiel wird die Auflösung (=Frequenz) mit 30 Minuten definiert durch "30T" und die Mittelwerte in diesem Zeitintervall berechnet. Die Formatierung ist definiert durch DateOffset. Anstatt des Mittelwerts könnte auch der Minimalwert, Maximalwert, Summenwert, u.s.w. berechnet werden - resample ist ein sehr mächtiges Instrument.

Zum Rechnen verwenden wir Numpy-Array welche über values vom DataFrame auslesen werden.

In [11]:
## Daten importieren
df = pd.read_csv('2017_DataExport15min.csv', delimiter=';', decimal='.', header = 0)
df['Time'] = pd.to_datetime(df['Time'], format='%Y-%m-%d %H:%M:%S' ) 
df = df.set_index(['Time']) # Datetime als Index definieren
df = df.resample('30T').mean() # Die Daten werden "resampled" zu 30 min Werte

hGlo  = df['hGlo'].values    # [W/m2] 
hDif  = df['hDif'].values    # [W/m2]
tAmb  = df['Tamb'].values    # [°C]
pVer  = df['Pload'].values   # [W] Verbrauchsprofil
zapf  = df['Zapfung'].values # [l/15min] Profil für Warmwasserbezug

Für die Sonnenstandsberechnung benötigen wir einen Array mit laufenden Stunden "lfStd", beginnend mit Stunde 0 bei Jahresanfang um die Zeitreihenberechnung umzusetzten.

Nun kann es sein, dass die einglesenen Daten zu einem beliebigem Zeitpunkt beginnen, die laufenden Stunden sind aber mit 1.Januar definiert, deshalb berechnen wir die Zeitdifferenz (zDif) zum 1.Januar des Jahres vom ersten Messwert.

Die Zeitdifferenz (zDif) ist im Format Timedelta, welches wir mit "total_seconds)" in ein Array mit Zahlen umwandeln. Die Zahlen sind Sekunden, welche wir in Stunden umrechnen, jedoch ist es nach wievor ein pandas-Array und wir benötigen ein Numpy-Array. Die erfolgt mit "lfStd.values" oder alternativ mit "np.array(lfStd).

In [12]:
zDif = df.index - dt.datetime(df.index[0].year, 1, 1, 0)
lfStd = zDif.total_seconds() # Pandas-Array mit Sekunden
lfStd = lfStd/3600 # Pandas-Array mit Stunden
lfStd = lfStd.values # Numpy-Array mit Zahlenwerte. Alternativ: lfStd = np.array(lfStd)

deltaT = lfStd[1] - lfStd[0] # [h] Auflösung
tutc = df.index

Wir haben folgendes gelernt:

  1. Importieren von Daten mit Pandas inklusiv Datum und Zeit
  2. Umwandeln des Datum-Zeit-Text in ein maschinenlesbares Format
  3. Resample der Daten

Simulation von einem Speicher

Für die Simulation einer Batterie ist neben der Produktion durch PV auch der Verbrauch relevant. Deshalb ist bei den Messdaten ein Lastprofil (pLoad) aufgeführt.

Bei der Betrachtung von Produktion oder Verbrauch konnten die Zeitreihen als Array vektoriell berechnet werden, was die Simulation sehr effizient machte. Bei der Berechnung des Ladezustands einer Batterie ist der vorgängige Zustand relevant, d.h. es kann nicht mehr vektoriell gerechnet werden, sondern über eine Schleife.

Bei der Berechnung wird jeder einzelne Zeitschritt berechnet:

1. Berechnung Leistungsbedarf (Verbrauch abzüglich Produktion)

Prüfen ob Leistungsbedarf durch die Batterie zur Verfügung gestellt werden kann:

a.) Wenn ja wird die Batterie entladen

b.) Wenn nicht wird die Leistung aus dem Netz bezogen, abzüglich dem Rest in der Batterie

2. Berechnung Überschuss (Produktion abzüglich Verbrauch)

Prüfen ob Batterie nicht voll ist,

a.) Wenn nicht kann die Batterie geladen werden bis sie voll ist, der Rest wird ins Netz rückgespeist

b.) Wenn die Batterie voll ist, wird der Überschuss ins Netz rückgespeist.

Mit diesem Vorgehen wird die Netzbezugsenergie (wBez) und die Netzeinspeisung (wEin) berechnet (Energie = Leistung mal Zeit).

Spätestens hier ist es sehr wichtig, wenn als Kommentar die Einheit angegeben wird, um Klarheit zu haben ob Leistung [W] oder Energie [kWh] oder [Wh] gemeint ist!

In [13]:
### Parameterdefinition
PR          = 0.75 #[1]
pPVstc      = 6  #[kW]
batCap      = 6  #[kWh]

pPV      = hGlo*PR*pPVstc # [W] PV-Produktionsleistung

### Batterie Berechnung
wBez  = 0    # [kWh] Energie Netzbezug
wEin  = 0    # [kWh] Energie Netzeinspeisung   
wBat  = 0    # [kWh] Energie Batteriestand. Annahme das Batterie zu Beginn leer ist

pBedarf = pVer-pPV  # [W] Bedarf ist Verbrauch agzüglich Produktion, jedoch....
pBedarf[pBedarf<0]=0 #...darf der Wert nicht <0 sein (z.B. wenn PV grösser ist als Verbrauch)
for p in range(pPV.size):    
    if wBat >= pBedarf[p]*deltaT/1000:
        wBat = wBat - pBedarf[p]*deltaT/1000    # [kWh]
    else:
        wBez = wBez + pBedarf[p]*deltaT/1000 - wBat   # [kWh]
        wBat  = 0
    pUebersch = np.maximum(pPV[p]-pVer[p], 0) # [W]      
    if wBat <= batCap:
        eNichtladbar = np.maximum(pUebersch*deltaT/1000 - (batCap-wBat),0) #[kWh]
        wBat = wBat + pUebersch*deltaT/1000 - eNichtladbar                 #[kWh] 
        wEin = wEin + eNichtladbar
    else:
        wEin = wEin + pUebersch*deltaT/1000   # [kWh]

wVer = np.sum(pVer)*deltaT/1000 # [kWh] Verbrauch Jahresenergie
wPro = np.sum(pPV)*deltaT/1000  # [kWh] PV-Produktion Jahresenergie
rEV = (wPro-wEin)/wPro *100 # [%] Eigenverbrauchsquote
rA = (wVer-wBez)/wVer *100 # [%] Autarkiegrad
print('Eigenverbrauchsgrad ist {0:1.0f} %'.format(rEV))
print('Autarkiegrad ist {0:1.0f} %'.format(rA))
Eigenverbrauchsgrad ist 52 %
Autarkiegrad ist 71 %

Mit dem obigen Code-Beispiel haben wir die Jahresenergie berechnet. Möchten wir einen zeitlichen Verlauf des Ladens und Entladens darstellen an drei Tagen im Jahr, so wird der Code erweitert, markiert mit (***):

In [14]:
### es wird nur ein Ausschnitt betrachtet, nicht das gesamte Jahr
a  = int(105*24/deltaT)
e  = int(a+3*24/deltaT)
pPV      = hGlo[a:e]*PR*pPVstc  # [W]
pVer     = pVer[a:e]
tutcTeil = tutc[a:e]

# Initialisierung Array für Darstellung (***)
batStand      = np.zeros(tutcTeil.size) # Array für Batteriestand (***)
netzbezug     = np.zeros(tutcTeil.size) # Array für Netzleistungsbezug (***)
einspeisung   = np.zeros(tutcTeil.size) # Array für Einspeiseleistung (***)

### Batterie Berechnung
wBez  = 0    # [kWh] Energie Netzbezug
wEin  = 0    # [kWh] Energie Netzeinspeisung   
wBat  = 0    # [kWh] Energie Batteriestand. Annahme das Batterie zu Beginn leer ist

pBedarf = pVer-pPV  # [W] Bedarf ist Verbrauch agzüglich Produktion, jedoch....
pBedarf[pBedarf<0]=0 #...darf der Wert nicht <0 sein (z.B. wenn PV grösser ist als Verbrauch)
for p in range(pPV.size):
    if wBat >= pBedarf[p]*deltaT/1000:
        wBat = wBat - pBedarf[p]*deltaT/1000    # [kWh]
        pBezug = 0                           # [kW] (***)
    else:
        wBez = wBez + pBedarf[p]*deltaT/1000 - wBat   # [kWh]
        wBat  = 0
        pBezug = pBedarf[p] - wBat*1000/deltaT   # [W] (***)
    pUebersch = np.maximum(pPV[p]-pVer[p], 0) # [W]      
    if wBat <= batCap:
        eNichtladbar = np.maximum(pUebersch*deltaT/1000 - (batCap-wBat),0) #[kWh]
        wBat = wBat + pUebersch*deltaT/1000 - eNichtladbar                 #[kWh] 
        wEin = wEin + eNichtladbar
        pEinspeisung = eNichtladbar/deltaT*1000  # [W] (***)
    else:
        wEin = wEin + pUebersch*deltaT/1000   # [kWh]
        pEinspeisung = pUebersch              # [W] (***)
    
    # Berechnete Werte ans Array übergeben (***)
    batStand[p]  = wBat*1000      # (***)
    netzbezug[p]   = pBezug       # [W] (***)
    einspeisung[p] = pEinspeisung # [W] (***)
 
### Darstellung (***)
fig = plt.figure(6, figsize=(10,4))
plt.plot_date(tutcTeil, pVer,'k', label='Verbrauch [W]')
plt.plot_date(tutcTeil, pPV,'g', label='Produktion [W]')
plt.plot_date(tutcTeil, batStand,'b', label='Batteriestand [Wh]')
plt.plot_date(tutcTeil, netzbezug,'r:', label='Netzbezug [W]')
plt.plot_date(tutcTeil, einspeisung,'g:', label='Rückspeisung [W]')
plt.legend(loc='best')
plt.xlabel('laufende Stunden [h]')
plt.ylabel('Leistung [W]/Batteriestand [Wh]')
plt.show()

Thermischer Speicher

Mit der selben Logik kann ein thermischer Speicher berechnet werden. Es wird von einer vollständigen Durchmischung im Speicher ausgegangen, d.h. es gibt keine Schichtung. Bei einem Energieeintrag erhöht sich die Temperatur im Speicher, was in der Praxis einer Stufenladung entspricht und keiner Schichtladung, da von keiner Schichtung ausgegangen wird.

Der Energieinhalt berechnet sich über die Speichertemperatur tsp, Kaltwassertemperatur tkw, der Dichte $\rho$ und der spezifischen Wärmekapazität von Wasser cp = 4.182 kJ/(kg K):

$W = vol \cdot \rho \cdot cp \cdot (tsp - tkw)$

Der Speicher hat folgendes Steuerverhalten: bei einer minimalen Speichertemperatur von tMin = 40 °C, wird der Heizstab über das Netz versorgt, sodass ausreichend Warmwasser erzeugt wird. Ist die Nachladung über das Netz aktiviert so bleibt diese aktiv bis der Speicher eine Temperatur von tMaxNetz = 45 °C hat (Hysteresis). Wird über PV geladen so ist die Ausschaltschwelle höher bei tMax = 85 °C. Es kann so mehr Energie eingelagert werden. Zu berücksichtigen sind Speicherverluste von 23 W.

In [15]:
### Thermische Speicher Berechnung
vol      = 300 # [liter] Speicherinhalt
tkw      = 12  # [°C] Kaltwassertemperatur
tMin     = 40  # [°C] min. Speichertemp. Heizung ein
tMax     = 85  # [°C] max. Speichertemp. bei PV-Ladung
tMaxNetz = 45  # [°C] max. Speichertemp. bei Netz-Ladung
thermalLosses = 23 # [W] Speicherverluste 
pHS           = 2  # [kW] Leistung Heizstab

tSp  = 45 # [°C] Sim.Starttemperatur Speicher       
wBez = 0  # [kWh] Energie Netzbezug
netzLaden = False # Flag welches Netzladung anzeigt
tSpArray = np.zeros(zapf.size)   # [°C] Array für Speichertemperatur
pThermArray = np.zeros(zapf.size)# [W] Array Leistung Heizstab
pzapfArray = np.zeros(zapf.size) # [W] Array Zapf-Leistung     

for p in range(pPV.size):
    if tSp < tMax:
        pTherm = pPV[p] # [W]
    else:
        pTherm = 0 # [W]

    if tSp > tMaxNetz:
        netzLaden = False

    if tSp < tMin:
        netzLaden = True

    if netzLaden == True:
        wBez = wBez + pHS*deltaT # [kWh]
        pTherm = pHS*1000 + pTherm # [W]

    wSp   = vol * 4.182/3.6 * (tSp-tkw) # [Wh] Energieinhalt Speicher vorher
    pzapf = zapf[p]/deltaT * 4.182/3.6 * (tSp-tkw) # [W_therm] gezapfte Leistung
    wSp   = wSp + (pTherm - pzapf - thermalLosses)*deltaT  # [Wh] Energieinhalt Speicher nachher
    tSp = wSp/(vol*4.182/3.6) + tkw # [°C]  neue Speichertemperatur
    tSpArray[p] = tSp # [°C]
    pThermArray[p] = pTherm # [W]
    pzapfArray[p] = pzapf # [W]

Dieses Tutorial wurde mit Jupyter Notebook erstellt, welches im Anaconda-Paket enthalten ist. Markus Markstaler. 2021

Dynamische Diagramme mit Bokeh

Für eine Visalisierung eines Strahlungsverlauf übers Jahr, wäre es nützlich wenn in die Daten "gezoomt" werden könnten, ohne der Eingrenzung über die Indexe [a:e]. Dies ist möglich mit dem Paket "bokeh". Es funktionert ähnlich wie Matplotlib. Es wird eine html-Datei erzeugt (mittels "output_file('filename.html')"), welche das Diagramm enthält und Werkzeuge zum Verschieben, Zommen, usw.

In [16]:
from bokeh.io import show  
from bokeh.plotting import figure 

x = np.linspace(0,7,100)
y = np.sin(x)

p1 = figure()
p1.line(x,y)

show(p1)

Nun zur Maximalvariante der Bokeh-Darstellung:

Die Übergabe der Daten erfolgt direkt als Dataframe (kann auch als Array erfolgen, wie bei Matplotlib). Im Dataframe ergänzen wir die Uhrzeit als String ("tutcStr") in einem, durch uns, definierbaren Format, zuvor wurde noch eine Stunde dazugezählt um von der UTC-Zeit auf die Mitteleuropäische Zeit zu kommen.

tt steht für Tooltip und heisst das die Werte der Kennlinie ausgegeben werden an der Position wo die Maus über der Kennlinie ist. Die Daten können die Messwerte sein, mit Zugriff "@hGlo" oder der Wert aus der gezeichneten Kennlinie mit $y", welcher interpoliert wird.

Bei der Erzeugung des Diagramm-Objekt kann die Grösse eingestellt werden und die Tooltips.

Eine sehr schöne Funktion ist p1.legend.click_policy="hide" dabei können die Kurven ausgeblendet werden, wenn die jeweilige Legende angeklickt wird.

In [17]:
from bokeh.io import output_file, output_notebook, show  # Bokeh is for making charts in a browser
from bokeh.plotting import figure 
from bokeh.embed import components


# zum Dateframe ergänzen wir die Uhrzeit als Text
tmez = tutc + pd.Timedelta('1 hours') # Ursprüngliche Zeit in UTC wird in MEZ umgewandelt (+1)
tmezStr = tmez.strftime("%d.%m.%Y %H:%M")

df['tmez'] = tmez 
df['tmezStr'] = tmezStr

output_file('myHTML.html') # das Diagramm wird in einem separaten html-File exportiert
output_notebook()   # das Diagramm wirm unten im JupyterNotebook dargestellt

tt = [
      ('Zeit MEZ:', '@tmezStr'),
      ('Globalstrahlung Diagramm:', '$y'),
      ('Globalstragkzbg Messwert:', '@hGlo'),
      ]

p1 = figure(title='Strahlungsverlauf', plot_width=800, plot_height=600, 
            tooltips = tt,  
            x_axis_type='datetime')
p1.line(x='tmez',y='hGlo', source=df, legend_label = 'Global', color="green")
p1.line(x='tmez',y='hDif', source=df, legend_label = 'Diffus', color="blue")
p1.xaxis.axis_label = 'Zeit [MEZ]'
p1.yaxis.axis_label = 'Strahlung [W/m²]'
p1.legend.location = 'top_right'
p1.legend.click_policy="hide" # Kurve ein/ausschaltbar
p1.toolbar.logo = None # deaktivieren des Bokeh-Logo
show(p1)
Loading BokehJS ...