0 Daumen
215 Aufrufe

Ich habe mir in Matlab einen Levenberg-Marquardt-Algorithmus geschrieben. Die Funktion funktioniert auch und je nachdem was ich optimieren will auch richtig.

Im Moment versuche ich gerade eine nichtlineare Funktion mit einer linearen Funktion zu approximieren. Genau gesagt, ich versuche mir mittels der Funktionen einen Materialparameter zu schätzen. Und genau hier liegt mein Problem. Irgendwie schafft es der Algorithmus nicht auf ein sinnvolles Ergebnis zu kommen. Ich bekomme immer kleinere Werte für meinen Parameter heraus.

Nun meine Frage: Kann ich denn eine untere Schranke in meinen Algorithmus einbauen, damit mein Parameter nicht kleiner werden kann als dieser Wert? Oder verfälscht das meine Optimierung?

Code:

function [a_est,v,error] = levmarq_unAd(predicted,a0,b0,x,un_1)

a0 = a0(:);
un_1 = un_1(:);

y_input = un_1;
v = zeros(1,30);
error = zeros(4,30);

Nparameters = length(a0);
n_iter = 30; % Anzahl der Iterationen
lambda = 0.01; % Anfangswert für den Dämpfungsfaktor
updateJ = 1;
h0 = 0.02; % Schrittweite

format long e;
a_est = a0;
a_est = a_est(:);

function r = res(x,a)
r = y_input - predicted(x,a);
r = r'*r;
end

for it = 1:n_iter

if updateJ == 1
% Auswertung der Jacobi-Matrix an den momentanen Parametern a_est
J = zeros(1, Nparameters);
E = eye(Nparameters);

% Auswertung des Totalen Fehlers (Abstand) an den momentanen
% Parametern
if it > 1
d = d_lm;
else
d = res(x,a_est);
end

for j = 1:Nparameters % spaltenweise durchlaufen
h = max(1, abs(a_est(j)) + 1e-2).*h0;
FDupper = res(x,(a_est + h.*E(:,j)));
FDlower = d;
J(:,j) = (1/h).*(FDupper - FDlower);
end

% Schrittweitenoptimierung
if norm(J(:,j))*h < norm(FDlower)*h0
h_neu = 0.5.*h;
J(:,j) = (1/h_neu).*(res(x,a_est + h_neu.*E(:,j)) - FDlower);
fprintf('h_neu = %.15f\n', h_neu);
end

% Berechnung der approximierten Hesse-Matrix, J' ist die
% Transponierte von J
H = J'*J;

if it == 1 % erste Iteration berechnet den Gesamtfehler
err = dot(d, d);
end

end

% fügt den Dämpfungsfaktor zur Hesse-Matrix hinzu
H_lm = H + (lambda.*diag(diag(H))); %eye(Nparameters, Nparameters));

% Berechnung der aktualisierten Parameter
dp = -inv(H_lm)*(J'*d(:));
a_lm = a_est + dp(:);

% Auswertung des Totalen Fehlers (Abstand) an den neuen Parametern
d_lm = res(x,a_lm);
err_lm = dot(d_lm, d_lm);

% Ist der Fehler der aktualisierten Parameter kleiner als der vorherige, dann
% mach sie zu den momentanen Parametern und reduziere den Dämpfungsfaktor
if err_lm < err
var = abs(a_lm - a_est);

if var < 1e-4
return;
end

if a_lm >= 500
lambda = lambda/10;
a_est = a_lm;
v(:,it) = a_est;
error(1,it) = err_lm;
error(2,it) = err;
error(3,it) = lambda;
error(4,it) = h;
err = err_lm;
fprintf('err = %.15f\n', err);
fprintf('a_est = %.15f\n', a_lm);
fprintf('it = %.15f\n', it);
updateJ = 1;
end

else % sonst erhöhe den Dämpfungsfaktor
updateJ = 0;
lambda = lambda.*10;
v(:,it) = a_est;
error(1,it) = err_lm;
error(2,it) = err;
error(3,it) = lambda;
error(4,it) = h;
fprintf('a_est = %.15f\n', a_lm);
fprintf('it = %.15f\n', it);
end

end

end

Ich hoffe das verwirrt jetzt nicht noch mehr. die if-Anweisung (a_lm >= 500) habe ich jetzt neu eingefügt, sozusagen als untere Schranke. Geht das? Oder muss ich das irgendwie anders lösen?


Nachtrag:

ich hab mir jetzt noch einmal ein paar gedanken gemacht... das mit der unteren Schranke ist natürlich Quatsch. Ich kann ja nicht bei einem lokalen minimierungsverfahren einfach werte "abschneiden"... mittlerweile vermute ich, dass es mit der fehlerberechnung zusammenhängt. also, dass ich wohl den fehler"falsch" berechne. im moment mache ich das ja über r = r'*r bzw. über err = dot(d, d);

Müsste ich eine andere Fehlernorm verwenden? Weiß da vielleicht jemand näheres?

von

Ein anderes Problem?

Stell deine Frage

Willkommen bei der Mathelounge! Stell deine Frage sofort und kostenfrei

x
Made by a lovely community