Ziel

Im Tutorial 2 von www.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.

Tutorial 2 - Ü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 erre icht die Strahlung einen Wert von 650 W/m2.

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

# Funktionsdefinition (benötigen wir später)
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)

Strahlungsarray (Zahlenreihe für die Strahlung in W/m2). Beachte das die Funktion sin() oben mit dem lambda-Operator definiert wurde und in Grad angegeben wird. Dies ist nicht bei der Funktion np.sin() der Fall.

Ü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 [29]:
np.shape(t)
Out[29]:
(1000,)

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

In [30]:
plt.figure(2, figsize=(8,4))
plt.plot(t,H, 'g')
plt.xlabel('Zeit [h]')
plt.ylabel('Strahlung [W/m2]')
Out[30]:
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".

In [31]:
A = 1 # Fläche 1 m2
deltaT = t[1] - t[0]
wH = np.sum(H)*deltaT*A/1000
print('Strahlungsenergie [kWh/m2/Tag]: ' + str(wH))
Strahlungsenergie [kWh/m2/Tag]: 4.965630132215672

Nun haben wir folgende Punkte gelernt:

  1. Erstellen eines Arrays.
  2. Darstellen von Arrays, inklusiv Beschriftung.
  3. Berechnen eines Flächenintegrals aus diskreten Werten, d.h. Aufsummieren von Werten.

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.

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

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

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

Tutorial 2 - Ü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.

Zum Importieren verwenden wir die Bibliothek pandas. Diese zeichnet sich durch gute Verarbeitung von Daten aus. Das Einlesen geht mittels "read_csv", dabei erhalten wir ein DataFrame (df). Dies ist vergleichbar mit einem Excel-Tabellenblatt mit einem Header, d.h. Überschrizeile und Links einem laufenden Index. Mittels "set_index" legen wir den Zeitstempel als Index fest.

Den Anfang desn DateFrames stellen wir mit "df.head()" dar.

Die csv-Datei ist auf github https://github.com/markstaler/pv4ing.

In [50]:
## Daten importieren
import pandas as pd

df = pd.read_csv('2017_DataExport15min.csv', delimiter=';', header = 0)
df = df.set_index(pd.DatetimeIndex(df['Time'])) 
print(df.head()) # zeigt den Header

hGlo  = df['hGlo'].values    # [W/m2] 
hDif  = df['hDif'].values    # [W/m2]
tAmb  = df['Tamb'].values    # [W/m2]
pVer  = df['Pload'].values   # [W] Verbrauchsprofil
zapf  = df['Zapfung'].values # [l/15min] Profil für Warmwasserbezug
tutc = df.index
                                    Time  hGlo  hDif  Tamb       Pload  \
Time                                                                     
2017-01-15 00:00:00  2017-01-15 00:00:00   0.0   0.0  -2.9  191.000000   
2017-01-15 00:15:00  2017-01-15 00:15:00   0.0   0.0  -3.0  100.000000   
2017-01-15 00:30:00  2017-01-15 00:30:00   0.0   0.0  -2.8  100.000000   
2017-01-15 00:45:00  2017-01-15 00:45:00   0.0   0.0  -3.0  100.000000   
2017-01-15 01:00:00  2017-01-15 01:00:00   0.0   0.0  -2.8  191.144444   

                         Wind  Niederschlag   hGloCOSMO   TambCOSMO  \
Time                                                                  
2017-01-15 00:00:00  0.640222           0.0         0.0   -2.531800   
2017-01-15 00:15:00  0.864444           0.0         0.0   -2.631250   
2017-01-15 00:30:00  0.495593           0.0         0.0   -2.730700   
2017-01-15 00:45:00  0.229046           0.0         0.0   -2.837075   
2017-01-15 01:00:00  0.375472           0.0         0.0   -2.943450   

                      NiederCOSMO  Zapfung  
Time                                        
2017-01-15 00:00:00           0.0      0.0  
2017-01-15 00:15:00           0.0      0.0  
2017-01-15 00:30:00           0.0      0.0  
2017-01-15 00:45:00           0.0      0.0  
2017-01-15 01:00:00           0.0      0.0  

Für die Sonnenstandsberechnung, später, erzeugen wir nun einen Array mit laufenden Stunden, beginnend mit Stunde 0 bei Jahresamfang.

In [51]:
lfStd = np.zeros(tutc.size)
for t in range(tutc.size):
    # berechnet laufender Tag im Jahr
    noDay = (tutc[t] - dt.datetime(tutc[0].year, 1, 1, 0)).days 
    # [h] berechnet laufende Stunde im Tag
    noHou = tutc[t].hour + (tutc[t].minute)/60.0 + (tutc[t].second)/3600.0 
    lfStd[t] = noDay*24 + noHou
    
deltaT = lfStd[1] - lfStd[0] # [h]

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 tutc[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 [52]:
hDir = hGlo - hDif

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

a=int(0)           # Anfangstag
e=int(7*24/deltaT) # Endtag

t = tutc[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[52]:
<matplotlib.legend.Legend at 0x236811eb630>

Dynamische Diagramme mit Bokeh

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

Zuerst noch die nützliche Funktion von pandas "resample", mit welcher eine Zeitreihe mit neuen Zeitschritte umgerechnet werden kann. Anschliessend werden die Daten mit "describe" dargestellt.

Die Übergabe der Daten erfolgt direkt als Dataframe (kann auch als Array erfolgen, wie bei Matplotlib). Im Dataframe erfänzen wir noch die Direktstrahlung und die Uhrzeit als String ("tutcStr") in einem, durch un, definierbaren Format.

Es werden die einzelnen Werkzeuge parametriert. "pan, box_zoom, reset, save" verstehen sich von selbst. Mit "hover" kann ein Tooltip ein/ausgeschalten werden. Tooltip heisst das die Daten 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 Kennlinie mit $y", welcher interpoliert wird.

Bei der Erzeugung des Diagramm-Objekt kann die Grösse eingestellt werden und die Parameter für Werkzeuge (tools) und Datenausgabe (tooltips).

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

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

## Daten importieren
df = pd.read_csv('2017_DataExport15min.csv', delimiter=';', header = 0)
df = df.set_index(pd.DatetimeIndex(df['Time'])) # erzeugt Datetime-Index welcher später für resampling notwendig ist
df = df.resample('30T').mean() # Die Daten werden "resampled" zu 30 min Werte. Dabei wird Mittelwertbildung durchgefürht
print(df.describe()) # gibt übersicht über die Zahlenwerte

tutc    = df.index # Zeit in der Datei ist UTC und wird in MEZ umgewandelt (+1)
tutcStr = tutc.strftime("%d.%m.%Y %H:%M")
    
df['tutc'] = tutc
df['tutcStr'] = tutcStr
df['hDir'] = df['hGlo'].values - df['hDif'].values

output_file('myHTML.html')
werkzeuge =  "pan, box_zoom, reset, save, hover"
tt = [
      ('Zeit UTC:', '@tutcStr'),
      ('Strahlung:', '$y'),
      ('Global:', '@hGlo'),
      ]

p1 = figure(title='Strahlungsverlauf', plot_width=1100, plot_height=600, 
            tooltips = tt, 
            tools = werkzeuge, 
            x_axis_type='datetime')
p1.line(x='tutc',y='hGlo', source=df, legend = 'Global', color="green")
p1.line(x='tutc',y='hDif', source=df, legend = 'Diffus', color="blue")
p1.line(x='tutc',y='hDir', source=df, legend = 'Direkt', color="red")
p1.xaxis.axis_label = 'Zeit'
p1.yaxis.axis_label = 'Strahlung [W/m²]'
p1.legend.location = 'top_right'
p1.legend.click_policy="hide" # Kurve ein/ausschaltbar
show(p1)
               hGlo          hDif          Tamb         Pload          Wind  \
count  17520.000000  17520.000000  17520.000000  17520.000000  17520.000000   
mean     137.970652     60.311662     11.307164    455.518114      1.076742   
std      230.948102     98.606569      8.593962    351.241692      1.005352   
min        0.000000      0.000000    -14.300000     77.950000      0.000000   
25%        0.000000      0.000000      4.650000    131.994444      0.539829   
50%        2.207416      0.039870     11.150000    468.572222      0.800349   
75%      190.916273     89.242393     17.600000    705.572222      1.193871   
max     1747.125753    923.758077     36.600000   1881.705555     30.863850   

       Niederschlag     hGloCOSMO     TambCOSMO   NiederCOSMO       Zapfung  
count  17520.000000  17520.000000  17520.000000  17520.000000  17520.000000  
mean       0.085537    125.915288     11.339318      0.143589      2.081570  
std        0.717983    203.829980      8.982928      0.538173      4.930293  
min        0.000000      0.000000    -15.299250      0.000000      0.000000  
25%        0.000000      0.000000      4.560537      0.000000      0.000000  
50%        0.000000      2.121812     10.987875      0.000000      0.000000  
75%        0.000853    186.890625     17.862594      0.006866      2.872567  
max       59.420356    873.362500     36.146000     10.418513     53.579618  

Schritt 3

Nun berechnen wir die Strahlungsenergie über das gesamte Jahr.

In [56]:
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('Jahrl. Energie Globalstrahlung ' + str(wGlo) + ' kWh/m2 entspricht ' + str(wGlo/wGlo*100) + ' %')
print('Jahrl. Energie Direktstrahlung ' + str(wDir) + ' kWh/m2 entspricht ' + str(wDir/wGlo*100) + ' %')
print('Jahrl. Energie Diffusstrahlung ' + str(wDif) + ' kWh/m2 entspricht ' + str(wDif/wGlo*100) + ' %')
Jahrl. Energie Globalstrahlung 1208.6229142200525 kWh/m2 entspricht 100.0 %
Jahrl. Energie Direktstrahlung 687.0434813297603 kWh/m2 entspricht 56.84514775008404 %
Jahrl. Energie Diffusstrahlung 528.3301597765812 kWh/m2 entspricht 43.71339923813399 %

Zur Kontrolle vergleiche die Jahresglobalstrahlungsenergie für diesen Standort mit den Werten aus PVGIS. Um anstatt der Jahresstrahlungsenergie die elektrische Energie aus einem Modul zu berechnen ist der Modulwirkungsgrad und die Modulfläche notwendig. Es muss noch der Temperatureinfluss und die Reflexionsverluste berücksichtigt werden. Theorie hierfür ist im Manuskript "Photovoltaik für Ingenieure" Kapitel "Kennwerte Photovoltaik".

Nun haben wir folgende Punkte gelernt:

  1. Bereinigen eines Array (z.B. Nullsetzten bei negativen Werten).
  2. Importieren von Messdaten.
  3. Berechnen der Strahlungskomponente hDir aus hGlob und hDiff.
  4. Darstellen der Kurven in einem Plot-Figure über eine Zeitachse, inklusiv einer Legende.
  5. Berechnen der JahresENERGIE aus der Jahresleistung über deltaT und aufsummieren.

Tutorial 2- Ü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

Hierfür wird das Dreikomponentenmodell verwendet, welches den Sonnenstand benötigt. Wie das Dreikomponentenmodell und der Sonnenstand berechnet wird, ist im Manuskript "Photvoltaik für Ingenieure" beschrieben. Hier der Code zur Berechnung der Jahresstrahlungsenergie bei unterschiedlicher Modul-Neigung und Azimut = 0° (Süden). 7 Code-Zeilen für den Sonnenstand und 7 Code-Zeilen für das Strahlentransformationsmodell.

In [57]:
# Parameterdefinition
lg = 9.4     # [grad] Längengrad Buchs
bg = 47.17   # [grad] Breitengrad Buchs
al = 0.2     # [1] Albedo-Faktor
aziFl = 0    # [grad] Azimut Modulfläche

# Sonnenstand
dekl = 23.45*cos(360/(365*24)*(lfStd - 173*24))     # [grad] Deklination
omega = 15 * lfStd + lg - 180                       # [grad] Stundenwinkel laufend
omega = np.mod(omega + 180, 360) - 180              # [grad] Stundenwinkel -180° bis 180°
h = asin(sin(dekl)*sin(bg) + cos(dekl)*cos(bg)*cos(omega)) # [grad] Sonnenhöhe
azi = acos((sin(h)*sin(bg) - sin(dekl))/(cos(h)*cos(bg)))
azi = np.sign(omega)*azi

h[h<0] = 0
neigung = np.linspace(0,90,50) # [grad] Neigung
wH = np.zeros(neigung.size) # Initialisierung
for n in range(neigung.size):
    neig = neigung[n]
    
    # Strahlentransformationsmodell
    cosTheta = cos(neig)*sin(h) + sin(neig)*cos(h)*cos(azi - aziFl)
    cosTheta[cosTheta<0] = 0
    hDifFl = hDif*(1 + cos(neig))/2
    hAlbFl = hGlo*al*(1 - cos(neig))/2
    hDirFl = hDir/np.maximum(sin(h), sin(5))*cosTheta
    hFl = hDirFl + hDifFl + hAlbFl
    wH[n] = np.sum(hFl)*deltaT/1000 # [kWh/m2] Jahresstrahlungsenergie    
    
# Darstellung
fig  = plt.figure(1, figsize=(8,5))
plt.plot(neigung, wH)
plt.xlabel('Neigung [°]')
plt.ylabel('Jahresstrahlungsenergie [°]')
plt.grid(which='both', linestyle='--')

Aufbau des Code-Beispiels: Für die unterschiedlichen Neigungen wird eine for-Schleife, um den Berechnungscode fürs Strahlentransformationsmodell gelegt. Innerhalb der for-Schleife wird die Leistung übers gesamte Jahr berechnet und für die Jahresenergie (wH) mit deltaT multipliziert und aufsummieren. Weiter wird die Summe durch 1000 geteilt um kWh anstatt Wh zu erhalten.

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

Zur Kontrolle:

  • Stelle den Sonnenstand dar (Azimut und Sonnenhöhe)
  • Vergleiche die Ergebnisse mit PVGIS für den berechneten Standort Buchs. 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: Variieren wir nicht einen Parameter, wie oben die Neigung, sondern einen zweiten Parameter, den Azimut, so lässt sich das Ergebnis als 3D-Plot darstellen oder als Contour-Plot. Der Code für ein 3D-Plot und Contour-Plot ist im Tutorial 1 beschrieben.

Zusammenfassung

Es wird für jeden Zeitpunkt eine Leistungsberechnung durchgeführt und am Ende des Jahres aufsummiert und mit dem Zeitschritt multipliziert. Somit wird die Jahresenergie berechnet (Energie = Leistung mal Zeit).

Die Berechnung der Jahresenergie erfolgt vektoriell.

Es werden nun mehrere Jahressimulationen durchgeführt für unterschiedliche Neigungen und für unterschiedliche Ausrichtungen, deshalb ist die Jahresberechnung in zwei for-Schleifen eingebettet.

Information zur 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. Das Lastprofil ist syntetisch erstellt aus folgenden Verbrauchern.

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 [58]:
### Parameterdefinition
PR          = 0.75 #[1]
pPVstc      = 2  #[kW]
batCap      = 2  #[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 '+str(rEV)+' %')
print('Autarkiegrad ist '+str(rA)+' %')
Eigenverbrauchsgrad ist 79.42460153785434 %
Autarkiegrad ist 36.08505487803508 %

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 [59]:
### 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 (Boiler) 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 [60]:
### 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 ebenfalls in Anaconda-Paket enthalten ist. Markus Markstaler.

In [ ]: