{
"cells": [
{
"cell_type": "markdown",
"id": "caa47946",
"metadata": {},
"source": [
"# 5. Funktionsfitting und Newtonverfahren (nD)\n",
"\n",
"Die PV-Anlage aus der vergangenen Übung soll nun als Modell abgebildet und als Modul eines übergeordneten Programms simuliert werden.\n",
"\n",
"Für die Modellbildung stehen bereits verschiedene Messwerte zur Verfügung sowie eine den Messwertverlauf beschreibende e-Funktion. Damit diese Funktion aber auch realistische Werte wiedergeben kann, müssen zuvor ihre Parameter optimiert werden. Dazu wird das mehrdimensionale Newtonverfahren angewendet.\n",
"\n",
"## Situation\n",
"\n",
"Zur Simulation einer PV-Anlage wurden normalisierte Leistungsdaten ermittelt. Die Messwerte wurden dazu über einen Zeitraum von 24h im 10min-Takt aufgezeichnet.\n",
"\n",
"Um das Modell der Anlage in ein übergeordnetes Programm einzubinden, muss aus den Messwerten nun noch eine Funktionsvorschrift abgeleitet werden (sog. Fitting).\n",
"\n",
"Dazu eignet sich folgende Funktion\n",
"\n",
"$$P(t) = e^{-\\left(\\frac{t-a}{b}\\right)^2}$$\n",
"\n",
"Zur Ermittlung der Funktionsparameter wird das mehrdimensionale Newtonverfahren verwendet und in MATLAB implementiert.\n",
"\n",
"PV_measure_A.txt\n",
"\n",
"PV_measure_B.txt\n",
"\n",
"## Aufgabe\n",
"\n",
"1. Importieren Sie Datensatz A und legen Sie $t$ und $P$ je als globale Variablen an."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "796f316d",
"metadata": {},
"outputs": [],
"source": [
"% your code here"
]
},
{
"cell_type": "markdown",
"id": "62365ea1",
"metadata": {},
"source": [
"2. Legen Sie eine neue Function `fit_fcn` an, in welcher die zu fittende Funktion abhängig von $t$ und den Parametern $a, b$ ausgewertet wird."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6e89e242",
"metadata": {},
"outputs": [],
"source": [
"function y = fit_fcn(t,params)\n",
" \t% “params” muss als Vektor übergeben werden\n",
" \ta = params(1);\n",
"\tb = params(2);\n",
" \ty = exp( -((t-a)/b)^2 );\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "8f273ccd",
"metadata": {},
"source": [
"3. Die Güte des Fittings wird durch den mittleren quadratischen Abstand zwischen Messwerten $P_m (t)$ und Funktionswerten $P(t)$ bestimmt. Die Funktionswerte hängen dabei von den gewählten Parametern $a, b$ ab:\n",
"\n",
"$$G(a,b) = \\frac{1}{n} \\sum_{i=1}^n {\\left(P(t_i)-P_{m_i}\\right)^2}$$\n",
"\n",
"Die Parameter $a, b$ sind also so zu bestimmen, dass $G$ minimal wird.\n",
"\n",
"Dazu muss nach beiden Parametern abgeleitet und Null gesetzt werden:\n",
"\n",
"$$\\frac{\\partial G}{\\partial a} = 0 = \\frac{1}{n}\\sum_{i=1}^n{2\\left(P(t_i)-P_{m_i}\\right)\\cdot\\frac{\\partial P(t_i)}{\\partial a}}$$\n",
"$$\\frac{\\partial G}{\\partial b} = 0 = \\frac{1}{n}\\sum_{i=1}^n{2\\left(P(t_i)-P_{m_i}\\right)\\cdot\\frac{\\partial P(t_i)}{\\partial b}}$$\n",
"\n",
"Zur Vereinfachung können die Konstanten $n$ und $2$ heraus gekürzt werden:\n",
"\n",
"$$0 = \\sum_{i=1}^n{\\left(P(t_i)-P_{m_i}\\right)\\cdot\\frac{\\partial P(t_i)}{\\partial a}}$$\n",
"$$0 = \\sum_{i=1}^n{\\left(P(t_i)-P_{m_i}\\right)\\cdot\\frac{\\partial P(t_i)}{\\partial b}}$$\n",
"\n",
"Sie benötigen also die partiellen Ableitungen der $e$-Funktion je nach $a$ und $b$. Diese Gleichungen werden **Bestimmungsgleichungen** genannt."
]
},
{
"cell_type": "markdown",
"id": "efe97d1d",
"metadata": {},
"source": [
"4. Legen Sie eine neue Function `fcn` an, in welcher Sie die Bestimmungsgleichungen abhängig von $a$ und $b$ auswerten:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "efa8cdc7",
"metadata": {},
"outputs": [],
"source": [
"function f = fcn(params)\n",
" \t% Aufruf der globalen Messwerte:\n",
" \tglobal tm Pm\n",
" \t% params ist ein Vektor und enthält a und b\n",
" \ta = params(1);\n",
"\tb = params(2);\n",
" \t% Vor der Summenbildung wird die e-Funktion mit den\n",
" \t% aktuellen Parametern ausgewertet\n",
" \tf0 = fit_fcn(tm,params);\n",
" \t% Nun werden die Bestimmungsgleichungen ausgerechnet:\n",
" \t% (Partielle Ableitungen einfügen!)\n",
" \tf(1,1) = sum( (f0-Pm).*%... );\n",
" \tf(2,1) = sum( (f0-Pm).*%... );\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "a7749ff2",
"metadata": {},
"source": [
"5. Um die Jacobi-Matrix einer Funktion numerisch zu ermitteln, können Sie folgende `function` verwenden ($f$ muss als Functionhandle `@(a,b)` übergeben werden):"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6036907b",
"metadata": {},
"outputs": [],
"source": [
"% Bestimmt die Jacobi-Matrix beliebiger Funktionen f durch Approximation\n",
"function J = jac(x,f)\n",
" % Dimension der Jacobi-Matrix bestimmt sich aus der Anzahl der Variablen \n",
" % (=Spalten) und Anzahl der Gleichungen (= Zeilen)\n",
" n = length(x);\n",
" % Funktionswerte an der Stelle x\n",
" f0 = f(x);\n",
" for i = 1:n\n",
" % Schrittweite zur Approximation der Ableitung\n",
" h = sqrt(eps)*max(1.0e-8,abs(x(i)));\n",
" % Zurücksetzen des x-Vektors auf Ausgangsposition\n",
" x1 = x;\n",
" % Erweitern der i-ten Komponenten von x um Schrittweite h\n",
" x1(i)=x1(i)+h;\n",
" % Berechnung einer Spalte der Jacobi-Matrix\n",
" J(:,i) = (f(x1)-f0)/h;\n",
" end\n",
"end"
]
},
{
"cell_type": "markdown",
"id": "5fc69e9e",
"metadata": {},
"source": [
"6. Die Vorschrift des mehrdimensionalen Newtonverfahrens lautet:\n",
"\n",
"$$J(x_k)\\cdot h_k = -F(x_k)$$\n",
"$$x_{k+1} = x_k + h_k$$\n",
"\n",
"Hierbei sind $F$ die Bestimmungsgleichungen und $J$ deren Jacobi-Matrix. Die Variable $x$ ist in unserem Fall ein Vektor und enthält die Parameter $a$ und $b$.\n",
"\n",
"Um das Newtonverfahren zu starten, benötigen Sie Startwerte für $a$ und $b$.\n",
"Es eignen sich z.B.:\n",
"\n",
"$a = 12.8, b = 3.5$\n",
"\n",
"Das lineare Gleichungssystem aus der ersten Zeile können Sie mit dem Backslash-Operator (`\\`) lösen (`mldivide`):"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "03442a2f",
"metadata": {},
"outputs": [],
"source": [
"A\\b = x;"
]
},
{
"cell_type": "markdown",
"id": "81834e86",
"metadata": {},
"source": [
"7. Werten Sie nach der Bestimmung der Parameter deren Güte aus. Wie groß ist der mittlere quadratische Abstand zwischen Messwerten und Fitting?"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "95e3a6a6",
"metadata": {},
"outputs": [],
"source": [
"% your code here"
]
},
{
"cell_type": "markdown",
"id": "ea3aaf6c",
"metadata": {},
"source": [
"8. Wenden Sie das Verfahren auch auf Datensatz B an.\n",
"\n",
"## Lösung\n",
"\n",
"\n",
"
- Hinweis: Dieses Video ist von YouTube aus eingebunden und nicht Teil des frei lizenzierten Materials! -
\n",
"\n",
"\n",
"Hier finden Sie außerdem den **Code** der Lösung:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ebeb02a3",
"metadata": {
"tags": [
"hide_input"
]
},
"outputs": [],
"source": [
"%% ---Hauptprogramm\n",
"% Im Programm wird das mehrdimensionale Newtonverfahren genutzt, um\n",
"% Leistungsdaten eine PV-Anlage möglichst gut durch eine\n",
"% Exponentialfunktion zu beschreiben.\n",
"function main\n",
" %--Impport der Messwerte\n",
" data = importdata('files/PV_measure_B.txt');\n",
" %--Anlegen globaler Variablen für Zeit (tm) und Leistung (Pm)\n",
" global tm Pm\n",
" tm = data(:,1);\n",
" Pm = data(:,2);\n",
" %--Plot der Messwerte\n",
" figure(1)\n",
" hold off\n",
" plot(tm,Pm)\n",
" grid on\n",
" %--Anlegen von Startwerten für das Newtonverfahren\n",
" % (Die Startwerte sind eine erste Schätzung der Parameter a und b)\n",
" params = [11;2.5];\n",
" %--Auswerten der zu fittenden Funktion mit den Startwerten über die\n",
" % gesamte Messdauer tm\n",
" P = fit_fcn(tm,params);\n",
" %--Berechnung der Norm der Bestimmungsgleichungen (soll gegen 0\n",
" % konvergieren)\n",
" norm(fcn(params))\n",
" %--Plotten des ersten Fitting-Versuchs\n",
" hold on\n",
" plot(tm,P)\n",
" %--Start des mehrdimensionalen Newtonverfahrens\n",
" % mit N Iterationsschritten\n",
" N = 200;\n",
" for i = 1:N\n",
" %--Die Variablen z0, aa und bb dienen nur der späteren Auswertung\n",
" % des Newtonverfahrens und sind nicht für die Lösung der Aufgabe\n",
" % notwendig\n",
" z0(i) = norm(fcn(params));\n",
" aa(i) = params(1);\n",
" bb(i) = params(2);\n",
" %--Berechnung der Ableitungswerte der Bestimmungsgleichungen durch\n",
" % Auswertung der Jacobimatrix\n",
" J = jac(params,@fcn);\n",
" %--Berechnung des Zwischenergebnisses h (Newton-\"Schrittweite\")\n",
" h = J\\-fcn(params);\n",
" %--Bestimmung eines neuen Parametersatz durch Addition von h auf\n",
" % den alten Parametersatz.\n",
" params = params + 0.5*h;\n",
" % Hinweis: Da das Newtonverfahren schwingen kann, wird es hier mit\n",
" % Faktor 0.5 gedämpft (0.5*h bzw. lambda*h mit lambda == 0.5).\n",
" % Das ist nicht immmer zwingend notwendig; Der Dämpfungsfaktor\n",
" % kann in anderen Fällen ggf. auch als lambda == 1 gewählt werden.\n",
" end\n",
" %--Finale Norm-Berechnung der Bestimmungsgleichungen (gibt Aufschluss\n",
" % über Güte des Ergebnis)\n",
" norm(fcn(params))\n",
" %--Bestimmung der Fitting-Werte über gesamten Messzeitraum\n",
" P = fit_fcn(tm,params);\n",
" %--Plotten der final gefitteten Exponentialfunktion\n",
" plot(tm,P,'k')\n",
" %--Beschriftung der Plots in Fenster 1\n",
" title('Funktionsfitting für eine PV-Anlage')\n",
" xlabel('Zeit [h]')\n",
" ylabel('Leistung (normalisiert)')\n",
" legend('Messwerte','Fitting 1 (Initial Guess)','Finaler Fit')\n",
" %--Ausgabe der finalen Funktionsparamter\n",
" params\n",
" %% Auswertung 1 des Newtonverfahrens\n",
" % (Diese Auswertung dient nur der Veranschaulichung des Verfahrens und\n",
" % ist nicht für die eigentliche Aufgabenstellung notwendig.\n",
" % Es werden die während dem Newtonverfahren gespeicherten\n",
" % Zwischenwerte der Parameter a und b, sowie die Norm der\n",
" % Bestimmungsgleichungen, geplottet.)\n",
" figure(2)\n",
" %--Norm der Bestimmungsgleichungen\n",
" subplot(2,1,1)\n",
" plot(z0)\n",
" grid on\n",
" title('Norm der Bestimmungsgleichungen')\n",
" xlabel('Iterationsschritte')\n",
" xlim([0 20])\n",
" %--Parameterwerte a/b\n",
" subplot(2,1,2)\n",
" plot(aa)\n",
" hold on\n",
" plot(bb)\n",
" hold off\n",
" grid on\n",
" title('Parameter a/b')\n",
" xlabel('Iterationsschritte')\n",
" xlim([0 20])\n",
" %% Auswertung 2 des Newtonverfahrens\n",
" % (Diese Auswertung dient nur der Veranschaulichung des Verfahrens und\n",
" % ist nicht für die eigentliche Aufgabenstellung notwendig.\n",
" % Sie sehen im entsprechenden Plot den Einfluss der Parameter a und b\n",
" % auf die Güte fes Funktionsfittings. Je näher der Z-Wert (=Norm der\n",
" % Bestimmungsgleichungen) gegen 0 geht, umso besser ist das Fitting.\n",
" % Das Newtonverfahren findet also ein Minimum in der über a und b\n",
" % aufgespannten Fläche.)\n",
" %--Auflösung bzw. Feinheit des Plots\n",
" N = N/4;\n",
" %--Definition von x/y-Gitterpunkten (für a und b)\n",
" xva = linspace(5,15,N);\n",
" xvb = linspace(0.5,5.5,N);\n",
" %--Berechung der jeweiligen z-Werte (=Norm der Bestimmungsgleichungen)\n",
" % abhängig von den \n",
" for i = 1:N\n",
" b = xvb(i);\n",
" for j = 1:N\n",
" a = xva(j);\n",
" Z(i,j) = norm(fcn([a,b]));\n",
" end\n",
" end\n",
" %--Anlegen der x/y-Werte (==a/b) als Matrix zum Plotten\n",
" [A,B] = meshgrid(xva,xvb);\n",
" %--Plotten der Normwerte über a und b als Surface-Plot\n",
" figure(3)\n",
" clf\n",
" hold off\n",
" surf(A,B,Z)\n",
" hold on\n",
" %--Plot des Verlaufs von z0 während dem Verfahren\n",
" plot3(aa,bb,z0,'r','LineWidth',2)\n",
"end\n",
"\n",
"%% ---Zu fittende Exponentialfunktion\n",
"% Funktion, welche die gemessenen Leistungsdaten der PV-Anlage möglichst\n",
"% gut wiedergeben soll\n",
"function y = fit_fcn(t,params)\n",
" a = params(1);\n",
" b = params(2);\n",
" y = exp( -((t-a)./b).^2 );\n",
"end\n",
"\n",
"%% ---Bestimmungsgleichungen für Funktionsparameter\n",
"% Sind notwendig, um mit dem Newotnverfahren die Werte für die Paramter der\n",
"% Exponentialfunktion zu optimieren (Güte wird durch mittleren quadrat.\n",
"% Abstand bestimmt)\n",
"function f = fcn(params)\n",
" % Aufruf globaler Variablen (Messwerte und Parameter p)\n",
" global tm Pm\n",
" % Parameter a,b,c und p werden im Vektor \"x\" zusammengefasst\n",
" a = params(1);\n",
" b = params(2);\n",
" % Werte der zu fittenden Funktion (hier: e-Funktion) an den Messwerten mit\n",
" % den entsprechenden Parametern (f(xi))\n",
" f0 = fit_fcn(tm,params);\n",
" % Bestimmungsgleichungen\n",
" f(1,1) = sum( (f0-Pm).*(2*(tm - a).*exp(-((tm - a)/b).^2))./b^2);\n",
" f(2,1) = sum( (f0-Pm).*(2*((tm - a).^2).*exp(-((tm - a)/b).^2))./b^3);\n",
"end\n",
"\n",
"%% ---Berechnung der Jacobi-Matrix\n",
"% Bestimmt die Jacobi-Matrix beliebiger Funktionen f durch Approximation\n",
"% Notwendig um die Ableitungswerte der Bestimmungsgleichungen zu erhalten\n",
"function J=jac(x,f)\n",
"% Dimension der Jacobimatrix bestimmt sich aus der Anzahl der Variablen \n",
"% (=Spalten) und Anzahl der Gleichungen (= Zeilen)\n",
"n = length(x);\n",
"% Funktionswerte an der Stelle x\n",
"f0 = f(x);\n",
"for i = 1:n\n",
" % Schrittweite zur Approximation der Ableitung\n",
" h = sqrt(eps)*max(1.0e-8,abs(x(i)));\n",
" % Zurücksetzen des x-Vektors auf Ausgangsposition\n",
" x1 = x;\n",
" % Erweitern der i-ten Komponenten von x um Schrittweite h\n",
" x1(i)=x1(i)+h;\n",
" % Berechnung einer Spalte der Jacobimatrix\n",
" J(:,i) = (f(x1)-f0)/h;\n",
"end\n",
"end"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Octave",
"language": "octave",
"name": "octave"
},
"language_info": {
"file_extension": ".m",
"help_links": [
{
"text": "GNU Octave",
"url": "https://www.gnu.org/software/octave/support.html"
},
{
"text": "Octave Kernel",
"url": "https://github.com/Calysto/octave_kernel"
},
{
"text": "MetaKernel Magics",
"url": "https://metakernel.readthedocs.io/en/latest/source/README.html"
}
],
"mimetype": "text/x-octave",
"name": "octave",
"version": "6.2.0"
}
},
"nbformat": 4,
"nbformat_minor": 5
}