Verfasst am: 21.12.2017, 08:47
Titel: Gauß-Newton-Algorithmus für nichtlineare Ausgleichsproblem
Hallo,
im Rahmen meines Studiums muss ich einen Gauß-Newton-Algorithmus zur Lösung von nichtlinearen Ausgleichsproblmen programmieren. Leider läuft mein bisher programmierter Quellcode noch nicht so ganz. Anbei ein Auszug der Aufgabenstellung:
Untersucht werden soll die nichtlineare Funktion:
𝑦 = 𝑓(𝑡) = 𝑎 + 𝑏 ⋅ 𝑒−𝑐𝑡
mit 3 Parametern a, b, c.
Die Parameter a, b, c sind so zu bestimmen, dass eine vorgegebene Punktwolke (Messwerte) (𝑡𝑖.𝑦𝑖), 𝑖 = 1, ,𝑛 nach dem Kriterium der Summe der kleinsten Fehlerquadrate durch die Funktion 𝑦 = 𝑓(𝑡) möglichst gut angenähert wird. Beispiele für Punktwolken mit n=10 und n=20 sind selbst zu erstellen. Programmierung des Gauß-Newton-Algorithmus für nichtlineare Ausgleichsprobleme mit matlab. Dabei soll die Anzahl der Punkte n variieren können. Der Startwert für die Parameter soll vorgegeben werden. Ein Abbruchkriterium ist zu definieren.
%Gauß-Newton-Algorthmus zur Lösung von nicht-linearen Ausgleichsproblemen
format long
tol = 1e-8; %Festlegung der Genauigkeit
maxstep = 30; %Festlegung des Abbruchkriteriums
t = [1,2,3,4,5,6,7,8,9,10]; %Generierung einer beliebigen Punktwolke
y = [2,4,6,8,2,3,5,7,4,2];
a = [1;1;1]; %Festlegung der Startwerte der Parameter a,b&c
m = length(t); %Bestimmung der Anzahl an Funktionpunkten
n = length(a); %Bestimmung der Anzahl von Parametern
aold = a;
for k=1:maxstep %Festlegung der Iteration
S = 0;
for i=1:m
for j=1:n %Anzahl der Parameter gibt die Anzahl an Spalten vor
J(i,j) = df(t(i),y(i),a(1,1),a(2,1),a(3,1),j) %Jacobimatrix
JT(j,i) = J(i,j) %Transponierte Jacobimatrix
end end
Jz = -JT*J; %Multiplikation der Jakobimatrix mit ihrer negiert-tranponierten
for i=1:m
r(i,1) = y(i)-a(1,1)+a(2,1)*exp(-a(3,1)*t(i)); %Aufstellen der Residuumsfunktion
S = S + r(i,1)^2 %Berechnung der Residuumsquadrate
end
S
g = Jz/JT;
a = aold-g*r; %Bewrechnung der neuen Approximation
unknowns = a;
err(k) = a(1,1)-aold(1,1); %Fehlerberechnung
if(abs(err(k)) <= tol); %Abbruch bei Unterschreitung der Tolleranz
break end
aold = a %aold entspricht der Lösung a nach definierter Genauigkeit
end
steps = k;
for i=1:m
f(i) = a(1,1)+a(2,1)*exp(-a(3,1)*x(i)); %Berechnung der Funktionswerte
end holdall plot(p,q, '*r'); %Plotten der Punktwolke
plot(a(1,1),a(2,1),a(3,1),'b'); %Plotten der Gauß-Newton-Approximation
title('Gauss-Newton-Approximation') %Titel,Legende und Achsenbeschriftung
xlabel('X') ylabel('Y') legend('Punktwolke','Gauss-Newton Approximation') hold off
end
function value = df(t,y,a,b,c,index) %Berechnung der partiellen Ableitung
switch index
case1
value = -b*c*exp(-c*t); %partielle Ableitung
end end
die fehlermeldung kommt daher, dass du deine funktion "df" mit dem index 2 aufrufst, da du aber nur case 1 definiert hast, wird der variable "value" kein wert zugewiesen.
ich würde mir den algorithmus nochmal genauer anschauen, auf wiki gibt es eine gute beschreibung, die man recht einfach nachvollziehen kann.
du hast meiner meinung einen fehler bei der berechnung der part. ableitungen, du bracuhst hier die ableitungen der residuumsfunktion nach den optimierungsparametern.
schonmal vielen Dank für deine schnelle Antwort. Habe mich noch mal über den Quellcode gesetzt und die Stellen korrigiert. Die Iteration läuft nun durch und ich bekomme mit dem derzeitigen Quellcode auch ein Ergebnis geplottet, jedoch scheint noch ein Fehler vorhanden zu sein-jedoch finde ich ihn einfach nicht....
%Gauß-Newton-Algorthmus zur Lösung von nicht-linearen Ausgleichsproblemen
format long
tol = 1e-8; %Festlegung der Genauigkeit
maxstep = 30; %Festlegung des Abbruchkriteriums
p = [1,2,3,4,5,6,7,8,9,10]; %Generierung einer beliebigen Punktwolke
q = [2,4,6,8,2,3,5,7,4,2];
a = [1;1.1;1.5]; %Festlegung der Startwerte der Parameter a,b&c
m = length(p); %Bestimmung der Anzahl an Funktionpunkten
n = length(a); %Bestimmung der Anzahl von Parametern
aold = a;
for k=1:maxstep %Festlegung der Iteration
S = 0;
for i=1:m
for j=1:n %Anzahl der Parameter gibt die Anzahl an Spalten vor
J(i,j) = df(p(i),a(1,1),a(2,1),a(3,1),j); %Jacobimatrix
JT(j,i) = J(i,j); %Transponierte Jacobimatrix
end end
Jz = -JT*J; %Multiplikation der Jakobimatrix mit ihrer negiert-tranponierten
for i=1:m
r(i,1) = (a(1,1)+a(2,1)*exp(-a(3,1)*p(i)))-q(i); %Aufstellen der Residuumsfunktion
S = r(i,1)^2 %Berechnung der Residuumsquadrate
end
g = Jz\JT;
a = aold-g*r; %Berechnung der neuen Approximation
unknowns = a;
err(k) = a(1,1)-aold(1,1); %Fehlerberechnung
if(abs(err(k)) <= tol); %Abbruch bei Unterschreitung der Tolleranz
break end
aold = a %aold entspricht der Lösung a nach definierter Genauigkeit
end
steps = k;
for i=1:m
f(i) = a(1,1)+a(2,1)*exp(-a(3,1)*p(i)) %Berechnung der Funktionswerte
end holdall plot(p,q, '*r'); %Plotten der Punktwolke
plot(p,f,'b'); %Plotten der Funktion bezüglich der Gauß-Newton-Approximation
title('Gauss-Newton-Approximation') %Titel,Legende und Achsenbeschriftung
xlabel('X') ylabel('Y') legend('Punktwolke','Gauss-Newton-Approximation') hold off
function value = df(p,a,b,c,index) %Berechnung der partiellen Ableitung
switch index
case1
value = -1;
case2
value = exp(-c*p);
case3
value = -p*b*exp(-c*p);
end
%Gauß-Newton-Algorithmus zur Lösung von nicht-linearen Ausgleichsproblemen
format long
tol = 1e-8; %Festlegung der Genauigkeit
maxstep = 30; %Festlegung des Abbruchkriteriums
p = [1,2,3,4,5,6,7,8,9,10]; %Generierung einer beliebigen Punktwolke
q = [2,2,2,8,2,3,5,7,4,2];
a = [2;5;7]; %Festlegung der Startwerte der Parameter a,b&c
m = length(p); %Bestimmung der Anzahl an Funktionpunkten
n = length(a); %Bestimmung der Anzahl von Parametern
aold = a;
for k=1:maxstep %Festlegung der Iteration
S = 0;
for i=1:m
for j=1:n %Anzahl der Parameter gibt die Anzahl an Spalten vor
J(i,j) = df(p(i),a(1,1),a(2,1),a(3,1),j); %Jacobimatrix
JT(j,i) = J(i,j); %Transponierte Jacobimatrix
end end
Jz = -JT*J; %Multiplikation der Jakobimatrix mit ihrer negiert-tranponierten
for i=1:m
r(i,1) = (a(1,1)+a(2,1)*exp(-a(3,1)*p(i)))-q(i); %Residuumsfunktion
S = S + r(i,1)^2 %Berechnung der Residuumsquadrate
end
g = Jz\JT;
a = aold-g*r; %Berechnung der neuen Approximation
unknowns = a;
err(k) = a(1,1)-aold(1,1); %Fehlerberechnung
if(abs(err(k)) <= tol); %Abbruch bei Unterschreitung der Tolleranz
break end
aold = a %aold entspricht der Lösung a nach definierter Genauigkeit
end
steps = k;
for i=1:m
f(i) = a(1,1)+a(2,1)*exp(-a(3,1)*p(i)) %Berechnung der Funktionswerte
end holdall plot(p,q, '*r'); %Plotten der Punktwolke
plot(p,f,'b'); %Plotten der Funktion bezüglich der Gauß-Newton-Approximation
title('Gauss-Newton-Approximation') %Titel,Legende und Achsenbeschriftung
xlabel('X') ylabel('Y') legend('Punktwolke','Gauss-Newton-Approximation') hold off
function value = df(p,a,b,c,index) %Berechnung der partiellen Ableitung
switch index
case1
value = 1;
case2
value = exp(-c*p);
case3
value = -p*b*exp(-c*p);
end
bin nun letztlich zu an einem, aus meiner Sicht passenden Ergebnis gelangt. Der Quellcode sieht dabei wie folgt aus:
Code:
% ------------------------------------------------------------------------ % Gauß-Newton-Algorthmus zur Lösung von nicht-linearen Ausgleichsproblemen % n = 20 Messpunkte % ------------------------------------------------------------------------ function[unknowns,steps,S] = GaussNewtonAlgorithmus() closeall% Schließung aller Fenster clearall% Leeren des Workspace clc% Leeren des Command Window format long
tol = 1e-4; % Toleranz
maxstep = 20; % Festlegung der max. Iteration
t = [1234567891011....% t-Werte der Punktwolke 121314151617181920];
y = [2221.51917.5161513.5... % y-Werte der Punktwolke 12.81211109.598.27.5...
7.16.76.56.25];
p = [2;21;0.5]; % Startwerte der Parameter a,b&c
m = length(t); % Anzahl an Zeilen
n = length(p); % Anzahl an Spalten
pold = p; % Erreichen der Genauigkeit
for k=1:maxstep % Festlegung der Iteration % ---------------------------------------------------------------------- % Bildung der Jacobi-Matrix und ihrer Transponierten % ---------------------------------------------------------------------- for i=1:m
for j=1:n
J(i,j) = df(t(i),p(1,1),p(2,1),p(3,1),j);
JT(j,i) = J(i,j);
end end % --------------------------------------------------------------------- % Aufstellen der Residuumsfunktion & Berechnung der Residuumsquadrate % ---------------------------------------------------------------------
S = 0;
for i=1:m
r(i,1) = (p(1,1)+p(2,1)*exp(-p(3,1)*t(i)))-y(i);
S = r(i,1)^2 end % --------------------------------------------------------------------- % Berechnung der Parameter über die Matrixgleichung % ---------------------------------------------------------------------
p = pold-(JT*J)\(JT*r);
unknowns = p;
% --------------------------------------------------------------------- % Erzeugung des Abbruchkriteriums bei Unterschreitung der Tolleranz % ---------------------------------------------------------------------
err(k) = p(1,1)-pold(1,1); % Fehlerberechnung if(abs(err(k)) <= tol);
break end
pold = p % Par. entsprechen Genauigkeit end
steps = k;
for i=1:m
f(i) = p(1,1)+p(2,1)*exp(-p(3,1)*t(i))% Berechnung der Funktionswerte end % ------------------------------------------------------------------------- % Grafische Darstellung der Ergebnisse % ------------------------------------------------------------------------- holdall;
grid on;
plot(t,y, '*r'); % Plotten der Punktwolke plot(t,f,'b'); % Plotten der Funktion % ------------------------------------------------------------------------- % Titel,Legende und Achsenbeschriftung % ------------------------------------------------------------------------- title('Gauss-Newton-Algorithmus') xlabel('t') ylabel('y') legend('data','best-fit') axis([021025]);
hold off
% ------------------------------------------------------------------------- % Ausgabe der Parameter a,b&c hinsichtlich "best fit" % -------------------------------------------------------------------------
a = p(1);
b = p(2);
c = p(3);
In = k;
str1 = ['Parameter a = ' num2str(p(1))];
str2 = ['Parameter b = ' num2str(p(2))];
str3 = ['Parameter c = ' num2str(p(3))];
str4 = ['Anzahl der Iterationen = ' num2str(k)];
disp(str1) disp(str2) disp(str3) disp(str4)
% ------------------------------------------------------------------------- % Switchbefehl für Parameter a,b&c / partielle Ableitungen % ------------------------------------------------------------------------- function value = df(t,a,b,c,index) switch index
case1
value = 1; % partielle Ableitung nach a case2
value = exp(-c*t); % partielle Ableitung nach b case3
value = -t*b*exp(-c*t); % partielle Ableitung nach c end
Du kannst Beiträge in dieses Forum schreiben. Du kannst auf Beiträge in diesem Forum antworten. Du kannst deine Beiträge in diesem Forum nicht bearbeiten. Du kannst deine Beiträge in diesem Forum nicht löschen. Du kannst an Umfragen in diesem Forum nicht mitmachen. Du kannst Dateien in diesem Forum posten Du kannst Dateien in diesem Forum herunterladen
MATLAB, Simulink, Stateflow, Handle Graphics, Real-Time Workshop, SimBiology, SimHydraulics, SimEvents, and xPC TargetBox are registered trademarks and The MathWorks, the L-shaped membrane logo, and Embedded MATLAB are trademarks of The MathWorks, Inc.
RSS