5. Funktionsfitting und Newtonverfahren (nD)

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

\[P(t) = e^{-\left(\frac{t-a}{b}\right)^2}\]

Zur Ermittlung der Funktionsparameter wird das mehrdimensionale Newtonverfahren verwendet und in MATLAB implementiert.

PV_measure_A.txt

PV_measure_B.txt

Aufgabe#

  1. Importieren Sie Datensatz A und legen Sie \(t\) und \(P\) je als globale Variablen an.

% your code here
  1. 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.

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
  1. 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:

\[G(a,b) = \frac{1}{n} \sum_{i=1}^n {\left(P(t_i)-P_{m_i}\right)^2}\]

Die Parameter \(a, b\) sind also so zu bestimmen, dass \(G\) minimal wird.

Dazu muss nach beiden Parametern abgeleitet und Null gesetzt werden:

\[\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}}\]
\[\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}}\]

Zur Vereinfachung können die Konstanten \(n\) und \(2\) heraus gekürzt werden:

\[0 = \sum_{i=1}^n{\left(P(t_i)-P_{m_i}\right)\cdot\frac{\partial P(t_i)}{\partial a}}\]
\[0 = \sum_{i=1}^n{\left(P(t_i)-P_{m_i}\right)\cdot\frac{\partial P(t_i)}{\partial b}}\]

Sie benötigen also die partiellen Ableitungen der \(e\)-Funktion je nach \(a\) und \(b\). Diese Gleichungen werden Bestimmungsgleichungen genannt.

  1. Legen Sie eine neue Function fcn an, 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
  1. Um die Jacobi-Matrix einer Funktion numerisch zu ermitteln, können Sie folgende function verwenden (\(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
  1. Die Vorschrift des mehrdimensionalen Newtonverfahrens lautet:

\[J(x_k)\cdot h_k = -F(x_k)\]
\[x_{k+1} = x_k + h_k\]

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;
  1. 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
  1. Wenden Sie das Verfahren auch auf Datensatz B an.

Lösung#

- Hinweis: Dieses Video ist von YouTube aus eingebunden und nicht Teil des frei lizenzierten Materials! -

Hier finden Sie außerdem den Code der Lösung:

Hide 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