Mein MATLAB Forum - goMatlab.de

Mein MATLAB Forum

 
Gast > Registrieren       Autologin?   
Bücher:

Fachkräfte:
Softwareentwickler Automotive Getriebesteuerungen (m/w)
Umsetzung der Softwarefunktionalität modellbasiert nach Spezifikation
MBtech Group GmbH & Co. KGaA - Sindelfingen

Systemingenieur (m/w) Funktionsentwicklung Automotive
Konzeption und Spezifikation von spezifischen Funktionen in elektronischen Steuergeräten
DRÄXLMAIER Group - Vilsbiburg bei Landshut

Applikationsingenieur (w/m) Testsysteme HIL-Simulation
Projektierung von Hardware-in-the-Loop-Systemen (Hardware und Software) in Kundenprojekten
dSPACE GmbH - Paderborn

Resident-Ingenieur (w/m) Hardware-in-the-Loop-Simulation
Inbetriebnahme und Software-Anpassungen der HIL-Systeme
dSPACE GmbH - Wolfsburg

System-Entwickler (m/w)
Analyse und Spezifikation von Systemanforderungen & Abstimmung von Anforderungen
MicroNova AG - Friedrichshafen

weitere Angebote

Partner:


Vermarktungspartner


Forum
      Option
[Erweitert]
  • Diese Seite per Mail weiterempfehlen
     


Gehe zu:  
Neues Thema eröffnen Neue Antwort erstellen

Gauß-Newton-Algorithmus für nichtlineare Ausgleichsprobleme

 

Schwalla
Forum-Newbie

Forum-Newbie


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


goMatlab ist ein Teil des goForen-Labels
goForen.de goMATLAB.de goLaTeX.de goPCB.de


 Impressum  | Nutzungsbedingungen  | Datenschutz  | Werbung/Mediadaten | Studentenversion | FAQ | goMatlab RSS Button RSS


Copyright © 2007 - 2018 goMatlab.de | Dies ist keine offizielle Website der Firma The Mathworks
Partner: LabVIEWforum.de

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.