3
Umweltsystemmodellierung (Fach) / Termin 4 Tutoratsaufgaben (Lektion)
- 1 siehe Abbildung Vergleich von Modell und Daten Die Datei „Zweispeicher.R“ enthält den Code für das Speichermodell mit zwei Speichern aus der Hausaufgabe der vergangenen Woche. Hier noch einmal die schematische Darstellung: Lesen Sie den Code der Funktion zweispeicher() mit dem Befehl source() ein (mit diesem Befehl können Sie R-Code aus einer Datei laden). Lesen Sie auch den Datensatz „Carnation_Creek_Daten.csv“ ein. Wenden Sie das Modell auf die Niederschlagsdaten vom 1. 8. 2007 bis 31. 7. 2008 an. Benutzen Sie die folgenden Parameterwerte: -Tabelle- Stellen Sie den zeitlichen Verlauf von gemessenem und modelliertem (Gesamt-)Abfluss grafisch dar. Versuchen Sie durch Ausprobieren die Parameterwerte so zu verändern, dass der modellierte Abfluss möglichst gut mit dem gemessenen übereinstimmt. # Funktion des Zwei-Speicher-Modells einlesen:source("Zweispeicher.R") # Carnation Creek - Daten einlesen:carnation <- read.csv("Carnation_Creek_Daten.csv") # Modell auf die Daten vom 1.8.2007 bis 31.7.2008 anwenden:tage <- as.Date(carnation$zeit, format = "%d.%m.%Y") begin <- which(tage=="2007-08-01")end <- which(tage=="2008-07-31")vorhersage <- zweispeicher(input=carnation$nied[begin:end], sa0=5, sb0=10, ka=0.9, kb=0.5, u=0.9) # Plot mit gemessenem und vorhergesagtem Abfluss:plot(tage[begin:end], carnation$qbeo[begin:end], xlab="Datum", ylab="Abfluss [mm/d]",type="l",col="black",las=1, lwd=2)lines(tage[begin:end], vorhersage$Abfluss_ges, col="red", lwd=2)legend("topright",legend=c("Daten","Modell"),lty=1,col=c("black","red"), lwd=2)
- 2. Kalibrierung a) Schreiben Sie eine Funktion, die den mittleren absoluten Fehler (MAF) des modellierten Gesamtabflusses im Vergleich zum gemessenen Abfluss berechnet. Als Argumente soll diese Funktion einen Vektor mit Niederschlagsdaten, einen Vektor mit gemessenen Abflussdaten und einen Vektor mit den Parameterwerten für die Funktion zweispeicher() bekommen. Zur Berechnung des modellierten Abflusses wird innerhalb der MAF-Funktion die zweispeicher-Funktion aufgerufen. Testen Sie die MAFFunktion mit den Daten und Parameterwerten aus Aufgabe 1. b) Benutzen Sie den Befehl optim(), um die Parameterwerte zu finden, für die der MAF mit den Daten vom 1. 8. 2007 bis 31. 7. 2008 („Trainingsdatensatz“) am geringsten ist. Benutzen Sie für die Optimierung die folgenden Startwerte: sa0=5, sb0=20, ka=0.8, kb=0.1, u=0.8. Lassen Sie sich die gefundenen optimalen Parameterwerte und den Wert des MAF für diese Parameterkombination ausgeben. c) Wenden Sie das zweispeicher-Modell mit der optimierten Parameterkombination auf den „Trainingsdatensatz“ an. Stellen Sie den zeitlichen Verlauf von gemessenem und modelliertem (Gesamt-)Abfluss grafisch dar. Vergleichen Sie die Grafik mit der aus Aufgabe 1. # a)# Funktion zur Berechnung des mittleren absoluten Fehlers# aus gemessenem Niederschlag und Abfluss und den Parametern der Modellfunktion: maf.speicher <- function(x, nied, abfl){model <- zweispeicher(input=nied, sa0=x[1], sb0=x[2], ka=x[3], kb=x[4], u=x[5])$Abfluss_gesfehler <- mean(abs(model-abfl))return(fehler)} # Test der Funktion:x <- c(5, 10, 0.9, 0.5, 0.9) # Parameterwertenied <- carnation$nied[begin:end]abfl <- carnation$qbeo[begin:end]maf.speicher(x, nied, abfl)# 4.591123 # b)# Parameteroptimierung:start <- c(5, 20, 0.8, 0.1, 0.8)optimal <- optim(par=start, fn=maf.speicher, nied=carnation$nied[begin:end], abfl=carnation$qbeo[begin:end]) # Ergebnis:optimal$par # die beste Parameterkombination# 3.356821e+00 2.203922e+01 4.930665e-01 5.474433e-04 5.890620e-01 optimal$value # der Wert der Zielfunktion für die beste Parameterkombination# 3.13334 # c)# Anwendung des Modells mit der besten Parameterkombination auf den "Trainingsdatensatz" (1.8.2007-31.7.2008):opt<-optimal$paroptvorhersage <- zweispeicher(input=carnation$nied[begin:end], sa0=opt[1], sb0=opt[2], ka=opt[3], kb=opt[4], u=opt[5]) plot(tage[begin:end], carnation$qbeo[begin:end], xlab="Datum", ylab="Abfluss [mm/d]",type="l",col="black",las=1, lwd=2)lines(tage[begin:end], optvorhersage$Abfluss_ges, col="red", lwd=2)legend("topright",legend=c("Daten","Modell"),lty=1,col=c("black","red"), lwd=2)
- 3. Evaluierung Wenden Sie das Modell mit den optimierten Parameterwerten (aus Aufgabe 2b) auf die Daten vom 1. 8. 2008 bis 31. 7. 2009 an. Stellen Sie auch hier den zeitlichen Verlauf von gemessenem und modelliertem (Gesamt-)Abfluss grafisch dar. Berechnen Sie den MAF für diesen Zeitraum. # Parameteroptimierung:start <- c(5, 20, 0.8, 0.1, 0.8)optimal <- optim(par=start, fn=maf.speicher, nied=carnation$nied[begin:end], abfl=carnation$qbeo[begin:end]) # Ergebnis:optimal$par # die beste Parameterkombination# 3.356821e+00 2.203922e+01 4.930665e-01 5.474433e-04 5.890620e-01 optimal$value # der Wert der Zielfunktion für die beste Parameterkombination# 3.13334 # c)# Anwendung des Modells mit der besten Parameterkombination auf den "Trainingsdatensatz" (1.8.2007-31.7.2008):opt<-optimal$paroptvorhersage <- zweispeicher(input=carnation$nied[begin:end], sa0=opt[1], sb0=opt[2], ka=opt[3], kb=opt[4], u=opt[5]) plot(tage[begin:end], carnation$qbeo[begin:end], xlab="Datum", ylab="Abfluss [mm/d]",type="l",col="black",las=1, lwd=2)lines(tage[begin:end], optvorhersage$Abfluss_ges, col="red", lwd=2)legend("topright",legend=c("Daten","Modell"),lty=1,col=c("black","red"), lwd=2) # Aufgabe 3: Evaluierung ----------------------------------------------------------------- # Grafische Darstellung von Modell und Daten für den Zeitraum 1.8.2008-31.7.2009:begin2 <- which(tage=="2008-08-01")end2 <- which(tage=="2009-07-31")opt<-optimal$paroptvorhersage2 <- zweispeicher(input=carnation$nied[begin2:end2], sa0=opt[1], sb0=opt[2], ka=opt[3], kb=opt[4], u=opt[5]) plot(tage[begin2:end2], carnation$qbeo[begin2:end2], xlab="Datum", ylab="Abfluss [mm/d]",type="l",col="black",las=1, lwd=2)lines(tage[begin2:end2], optvorhersage2$Abfluss_ges, col="red", lwd=2)legend("topright",legend=c("Daten","Modell"),lty=1,col=c("black","red"), lwd=2) # Berechnung des MAF für die neuen Daten:maf.speicher(opt, carnation$nied[begin2:end2], carnation$qbeo[begin2:end2])# 2.608325