WICHTIG: Der Betrieb von goMatlab.de wird privat finanziert fortgesetzt. - Mehr Infos...

Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   

Partner:




Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Gauß-Newton-Algorithmus für nichtlineare Ausgleichsproblem

 

Schwalla
Forum-Newbie

Forum-Newbie


Beiträge: 8
Anmeldedatum: 14.07.16
Wohnort: Hannover
Version: ---
     Beitrag Verfasst am: 21.12.2017, 08:47     Titel: Gauß-Newton-Algorithmus für nichtlineare Ausgleichsproblem
  Antworten mit Zitat      
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.

Code:
function [unknowns,steps,S] = GaussNewton()

%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
hold all
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
    case 1
    value = -b*c*exp(-c*t); %partielle Ableitung
end
end


Als Fehler wird mir im Command Window
"(line18)
J(i,j) = df(t(i),y(i),a(1,1),a(2,1),a(3,1),j) %Jacobimatrix" angezeigt.

LG
Private Nachricht senden Benutzer-Profile anzeigen


Friidayy
Forum-Century

Forum-Century


Beiträge: 225
Anmeldedatum: 17.12.13
Wohnort: ---
Version: R2012b
     Beitrag Verfasst am: 21.12.2017, 09:52     Titel:
  Antworten mit Zitat      
hallo schwalla,

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.

https://de.wikipedia.org/wiki/Gau%C3%9F-Newton-Verfahren

du hast meiner meinung einen fehler bei der berechnung der part. ableitungen, du bracuhst hier die ableitungen der residuumsfunktion nach den optimierungsparametern.

gruß, michael
Private Nachricht senden Benutzer-Profile anzeigen
 
Schwalla
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 8
Anmeldedatum: 14.07.16
Wohnort: Hannover
Version: ---
     Beitrag Verfasst am: 21.12.2017, 17:18     Titel:
  Antworten mit Zitat      
Hallo Michael,

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....
Code:
function [unknowns,steps,S] = GaussNewtonAlgorithmus()

%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
hold all
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
    case 1
    value = -1;
    case 2
    value = exp(-c*p);
    case 3
    value = -p*b*exp(-c*p);
end


Der Fehler liegt dabei laut matlab in line 30:

(line 30)
Warning: Matrix is singular, close to singular or badly scaled. Results may be inaccurate. RCOND =
NaN.

Confused Confused Confused Confused

LG
Private Nachricht senden Benutzer-Profile anzeigen
 
Friidayy
Forum-Century

Forum-Century


Beiträge: 225
Anmeldedatum: 17.12.13
Wohnort: ---
Version: R2012b
     Beitrag Verfasst am: 21.12.2017, 22:11     Titel:
  Antworten mit Zitat      
deine ansatzfunktion lautet: f(x) = a + b * exp (c * x), richtig?

ableitung nach a lautet dann 1 und nicht -1 (zeile 60). damit läuft das programm aber leider immer noch nicht rund.
Private Nachricht senden Benutzer-Profile anzeigen
 
Schwalla
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 8
Anmeldedatum: 14.07.16
Wohnort: Hannover
Version: ---
     Beitrag Verfasst am: 21.12.2017, 22:44     Titel:
  Antworten mit Zitat      
Hallo,

die Funktion lautet:

f(t) = a + b * exp (-c * t)

Aber du hast Recht, es muss natürlich nach a abgeleitet 1 sein.

Anbei nun der korrigierte immer noch fehlerbehaftete code Confused Confused Confused

Code:
function [unknowns,steps,S] = GaussNewtonAlgorithmus()

%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
hold all
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
    case 1
    value = 1;
    case 2
    value = exp(-c*p);
    case 3
    value = -p*b*exp(-c*p);
end


LG
Private Nachricht senden Benutzer-Profile anzeigen
 
Friidayy
Forum-Century

Forum-Century


Beiträge: 225
Anmeldedatum: 17.12.13
Wohnort: ---
Version: R2012b
     Beitrag Verfasst am: 01.01.2018, 14:54     Titel:
  Antworten mit Zitat      
lösung als referenz


Code:
function main
m.x = [1 2 3 4 5 6 7 8 9 10];
m.y = [2 4 6 8 2 3 5 7 4 2];
p0 = [0 0 0]';

p = fminsearch(@(p)minFunc(p,m),p0);

figure(1);
plot(m.x,m.y,'+'); hold on; plot(m.x,fitFunc(m.x,p));
end

function y = fitFunc(x,p)
y = zeros(size(x));
for i = 1:length(x)
    y(i) = p(1)+p(2)*exp(-p(3)*x(i));
end
end

function J = minFunc(p,m)
J = sum((m.y-fitFunc(m.x,p)).^2);
end
Private Nachricht senden Benutzer-Profile anzeigen
 
Schwalla
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 8
Anmeldedatum: 14.07.16
Wohnort: Hannover
Version: ---
     Beitrag Verfasst am: 03.01.2018, 10:14     Titel:
  Antworten mit Zitat      
Hallo,

vielen Dank hat mich schon mal weiter gebracht Wink
Private Nachricht senden Benutzer-Profile anzeigen
 
Schwalla
Themenstarter

Forum-Newbie

Forum-Newbie


Beiträge: 8
Anmeldedatum: 14.07.16
Wohnort: Hannover
Version: ---
     Beitrag Verfasst am: 20.01.2018, 20:17     Titel:
  Antworten mit Zitat      
Hallo,

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()
close all                                 % Schließung aller Fenster
clear all                                 % Leeren des Workspace
clc                                       % Leeren des Command Window
format long
tol = 1e-4;                               % Toleranz
maxstep = 20;                             % Festlegung der max. Iteration
t = [1 2 3 4 5 6 7 8 9 10 11....          % t-Werte der Punktwolke
     12 13 14 15 16 17 18 19 20];
y = [22 21.5 19 17.5 16 15 13.5...        % y-Werte der Punktwolke
     12.8 12 11 10 9.5 9 8.2 7.5...
     7.1 6.7 6.5 6.2 5];  
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
% -------------------------------------------------------------------------
hold all;
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([0 21 0 25]);
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
    case 1
    value = 1;                              % partielle Ableitung nach a
    case 2
    value = exp(-c*t);                      % partielle Ableitung nach b
    case 3
    value = -t*b*exp(-c*t);                 % partielle Ableitung nach c
end
 


Vielen Dank nochmals für die Hilfestellung Smile
LG
Private Nachricht senden Benutzer-Profile anzeigen
 
Neues Thema eröffnen Neue Antwort erstellen



Einstellungen und Berechtigungen
Beiträge der letzten Zeit anzeigen:

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
.





 Impressum  | Nutzungsbedingungen  | Datenschutz | FAQ | goMatlab RSS Button RSS

Hosted by:


Copyright © 2007 - 2024 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks

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.