5. Funktionsfitting und Newtonverfahren (nD)#
Die PV-Anlage aus der vergangenen Übung soll nun als Modell abgebildet und als Modul eines übergeordneten Programms simuliert werden.
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.
Situation#
Zur Simulation einer PV-Anlage wurden normalisierte Leistungsdaten ermittelt. Die Messwerte wurden dazu über einen Zeitraum von 24h im 10min-Takt aufgezeichnet.
Um das Modell der Anlage in ein übergeordnetes Programm einzubinden, muss aus den Messwerten nun noch eine Funktionsvorschrift abgeleitet werden (sog. Fitting).
Dazu eignet sich folgende Funktion
Zur Ermittlung der Funktionsparameter wird das mehrdimensionale Newtonverfahren verwendet und in MATLAB implementiert.
Aufgabe#
Importieren Sie Datensatz A und legen Sie \(t\) und \(P\) je als globale Variablen an.
% your code here
Legen Sie eine neue Function
fit_fcnan, in welcher die zu fittende Funktion abhängig von \(t\) und den Parametern \(a, b\) ausgewertet wird.
function y = fit_fcn(t,params)
% “params” muss als Vektor übergeben werden
a = params(1);
b = params(2);
y = exp( -((t-a)/b)^2 );
end
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:
Die Parameter \(a, b\) sind also so zu bestimmen, dass \(G\) minimal wird.
Dazu muss nach beiden Parametern abgeleitet und Null gesetzt werden:
Zur Vereinfachung können die Konstanten \(n\) und \(2\) heraus gekürzt werden:
Sie benötigen also die partiellen Ableitungen der \(e\)-Funktion je nach \(a\) und \(b\). Diese Gleichungen werden Bestimmungsgleichungen genannt.
Legen Sie eine neue Function
fcnan, in welcher Sie die Bestimmungsgleichungen abhängig von \(a\) und \(b\) auswerten:
function f = fcn(params)
% Aufruf der globalen Messwerte:
global tm Pm
% params ist ein Vektor und enthält a und b
a = params(1);
b = params(2);
% Vor der Summenbildung wird die e-Funktion mit den
% aktuellen Parametern ausgewertet
f0 = fit_fcn(tm,params);
% Nun werden die Bestimmungsgleichungen ausgerechnet:
% (Partielle Ableitungen einfügen!)
f(1,1) = sum( (f0-Pm).*%... );
f(2,1) = sum( (f0-Pm).*%... );
end
Um die Jacobi-Matrix einer Funktion numerisch zu ermitteln, können Sie folgende
functionverwenden (\(f\) muss als Functionhandle@(a,b)übergeben werden):
% Bestimmt die Jacobi-Matrix beliebiger Funktionen f durch Approximation
function J = jac(x,f)
% Dimension der Jacobi-Matrix bestimmt sich aus der Anzahl der Variablen
% (=Spalten) und Anzahl der Gleichungen (= Zeilen)
n = length(x);
% Funktionswerte an der Stelle x
f0 = f(x);
for i = 1:n
% Schrittweite zur Approximation der Ableitung
h = sqrt(eps)*max(1.0e-8,abs(x(i)));
% Zurücksetzen des x-Vektors auf Ausgangsposition
x1 = x;
% Erweitern der i-ten Komponenten von x um Schrittweite h
x1(i)=x1(i)+h;
% Berechnung einer Spalte der Jacobi-Matrix
J(:,i) = (f(x1)-f0)/h;
end
end
Die Vorschrift des mehrdimensionalen Newtonverfahrens lautet:
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\).
Um das Newtonverfahren zu starten, benötigen Sie Startwerte für \(a\) und \(b\). Es eignen sich z.B.:
\(a = 12.8, b = 3.5\)
Das lineare Gleichungssystem aus der ersten Zeile können Sie mit dem Backslash-Operator (\) lösen (mldivide):
A\b = x;
Werten Sie nach der Bestimmung der Parameter deren Güte aus. Wie groß ist der mittlere quadratische Abstand zwischen Messwerten und Fitting?
% your code here
Wenden Sie das Verfahren auch auf Datensatz B an.
Lösung#
Hier finden Sie außerdem den Code der Lösung:
Show code cell source
%% ---Hauptprogramm
% Im Programm wird das mehrdimensionale Newtonverfahren genutzt, um
% Leistungsdaten eine PV-Anlage möglichst gut durch eine
% Exponentialfunktion zu beschreiben.
function main
%--Impport der Messwerte
data = importdata('files/PV_measure_B.txt');
%--Anlegen globaler Variablen für Zeit (tm) und Leistung (Pm)
global tm Pm
tm = data(:,1);
Pm = data(:,2);
%--Plot der Messwerte
figure(1)
hold off
plot(tm,Pm)
grid on
%--Anlegen von Startwerten für das Newtonverfahren
% (Die Startwerte sind eine erste Schätzung der Parameter a und b)
params = [11;2.5];
%--Auswerten der zu fittenden Funktion mit den Startwerten über die
% gesamte Messdauer tm
P = fit_fcn(tm,params);
%--Berechnung der Norm der Bestimmungsgleichungen (soll gegen 0
% konvergieren)
norm(fcn(params))
%--Plotten des ersten Fitting-Versuchs
hold on
plot(tm,P)
%--Start des mehrdimensionalen Newtonverfahrens
% mit N Iterationsschritten
N = 200;
for i = 1:N
%--Die Variablen z0, aa und bb dienen nur der späteren Auswertung
% des Newtonverfahrens und sind nicht für die Lösung der Aufgabe
% notwendig
z0(i) = norm(fcn(params));
aa(i) = params(1);
bb(i) = params(2);
%--Berechnung der Ableitungswerte der Bestimmungsgleichungen durch
% Auswertung der Jacobimatrix
J = jac(params,@fcn);
%--Berechnung des Zwischenergebnisses h (Newton-"Schrittweite")
h = J\-fcn(params);
%--Bestimmung eines neuen Parametersatz durch Addition von h auf
% den alten Parametersatz.
params = params + 0.5*h;
% Hinweis: Da das Newtonverfahren schwingen kann, wird es hier mit
% Faktor 0.5 gedämpft (0.5*h bzw. lambda*h mit lambda == 0.5).
% Das ist nicht immmer zwingend notwendig; Der Dämpfungsfaktor
% kann in anderen Fällen ggf. auch als lambda == 1 gewählt werden.
end
%--Finale Norm-Berechnung der Bestimmungsgleichungen (gibt Aufschluss
% über Güte des Ergebnis)
norm(fcn(params))
%--Bestimmung der Fitting-Werte über gesamten Messzeitraum
P = fit_fcn(tm,params);
%--Plotten der final gefitteten Exponentialfunktion
plot(tm,P,'k')
%--Beschriftung der Plots in Fenster 1
title('Funktionsfitting für eine PV-Anlage')
xlabel('Zeit [h]')
ylabel('Leistung (normalisiert)')
legend('Messwerte','Fitting 1 (Initial Guess)','Finaler Fit')
%--Ausgabe der finalen Funktionsparamter
params
%% Auswertung 1 des Newtonverfahrens
% (Diese Auswertung dient nur der Veranschaulichung des Verfahrens und
% ist nicht für die eigentliche Aufgabenstellung notwendig.
% Es werden die während dem Newtonverfahren gespeicherten
% Zwischenwerte der Parameter a und b, sowie die Norm der
% Bestimmungsgleichungen, geplottet.)
figure(2)
%--Norm der Bestimmungsgleichungen
subplot(2,1,1)
plot(z0)
grid on
title('Norm der Bestimmungsgleichungen')
xlabel('Iterationsschritte')
xlim([0 20])
%--Parameterwerte a/b
subplot(2,1,2)
plot(aa)
hold on
plot(bb)
hold off
grid on
title('Parameter a/b')
xlabel('Iterationsschritte')
xlim([0 20])
%% Auswertung 2 des Newtonverfahrens
% (Diese Auswertung dient nur der Veranschaulichung des Verfahrens und
% ist nicht für die eigentliche Aufgabenstellung notwendig.
% Sie sehen im entsprechenden Plot den Einfluss der Parameter a und b
% auf die Güte fes Funktionsfittings. Je näher der Z-Wert (=Norm der
% Bestimmungsgleichungen) gegen 0 geht, umso besser ist das Fitting.
% Das Newtonverfahren findet also ein Minimum in der über a und b
% aufgespannten Fläche.)
%--Auflösung bzw. Feinheit des Plots
N = N/4;
%--Definition von x/y-Gitterpunkten (für a und b)
xva = linspace(5,15,N);
xvb = linspace(0.5,5.5,N);
%--Berechung der jeweiligen z-Werte (=Norm der Bestimmungsgleichungen)
% abhängig von den
for i = 1:N
b = xvb(i);
for j = 1:N
a = xva(j);
Z(i,j) = norm(fcn([a,b]));
end
end
%--Anlegen der x/y-Werte (==a/b) als Matrix zum Plotten
[A,B] = meshgrid(xva,xvb);
%--Plotten der Normwerte über a und b als Surface-Plot
figure(3)
clf
hold off
surf(A,B,Z)
hold on
%--Plot des Verlaufs von z0 während dem Verfahren
plot3(aa,bb,z0,'r','LineWidth',2)
end
%% ---Zu fittende Exponentialfunktion
% Funktion, welche die gemessenen Leistungsdaten der PV-Anlage möglichst
% gut wiedergeben soll
function y = fit_fcn(t,params)
a = params(1);
b = params(2);
y = exp( -((t-a)./b).^2 );
end
%% ---Bestimmungsgleichungen für Funktionsparameter
% Sind notwendig, um mit dem Newotnverfahren die Werte für die Paramter der
% Exponentialfunktion zu optimieren (Güte wird durch mittleren quadrat.
% Abstand bestimmt)
function f = fcn(params)
% Aufruf globaler Variablen (Messwerte und Parameter p)
global tm Pm
% Parameter a,b,c und p werden im Vektor "x" zusammengefasst
a = params(1);
b = params(2);
% Werte der zu fittenden Funktion (hier: e-Funktion) an den Messwerten mit
% den entsprechenden Parametern (f(xi))
f0 = fit_fcn(tm,params);
% Bestimmungsgleichungen
f(1,1) = sum( (f0-Pm).*(2*(tm - a).*exp(-((tm - a)/b).^2))./b^2);
f(2,1) = sum( (f0-Pm).*(2*((tm - a).^2).*exp(-((tm - a)/b).^2))./b^3);
end
%% ---Berechnung der Jacobi-Matrix
% Bestimmt die Jacobi-Matrix beliebiger Funktionen f durch Approximation
% Notwendig um die Ableitungswerte der Bestimmungsgleichungen zu erhalten
function J=jac(x,f)
% Dimension der Jacobimatrix bestimmt sich aus der Anzahl der Variablen
% (=Spalten) und Anzahl der Gleichungen (= Zeilen)
n = length(x);
% Funktionswerte an der Stelle x
f0 = f(x);
for i = 1:n
% Schrittweite zur Approximation der Ableitung
h = sqrt(eps)*max(1.0e-8,abs(x(i)));
% Zurücksetzen des x-Vektors auf Ausgangsposition
x1 = x;
% Erweitern der i-ten Komponenten von x um Schrittweite h
x1(i)=x1(i)+h;
% Berechnung einer Spalte der Jacobimatrix
J(:,i) = (f(x1)-f0)/h;
end
end