Verfasst am: 23.02.2013, 16:55
Titel: Optimierung/Minimierung mehrerer Parameter eines DGL-Systems
Hallo allerseits,
ich habe derzeit ein Problem und würde mich über einen Lösungsvorschlag und Unterstützung sehr freuen.
Ich habe ein Gleichungssystem aus 4 Differentialgleichungen mit 4 Unbekannten (x1-4) und 5 Parametern (a1-5).
Die Lösungen (oder Lösungsvektor, die ich mit Matlab errechne) stellen Modelldaten dar. ZUsätzlich habe ich auch reale Daten (Messdaten).
Nun sollen die Parametern a1-5 mit der "Methode der kleinsten Quadradrate" optimiert/geschätz werden, so das die Modelldaten
näherungsweise den Messdaten entsprechen.
Hier das DGL-System, was das System, aus dem ich Messdaten habe beschreibt:
dy1(t)/dt = -a1*x1(t)
dy2(t)/dt = -a2*x2(t) + a1*x1(t)
dy3(t)/dt = -a3*x3(t) + a2*X2(t) + a4*x4(t)
dy4(t)/dt = -a5*x4(t) + a3*x3(t) - a4*x4(t)
Mein Problem ist nun, das ich nicht weiss, wie ich das Problem mit Matlab/Octave (in der Uni arbeite ich mit MAtlab , Zuhause mit Octave, wenn es für Octave eine andere Lösung geben würde, wäre ich an der Octave Lösung interessiert) lösen kann, speziell, da es sich um ein DGL-System handelt.
Das Lösen des DGL-System ist kein Problem. Auch die Optimierung, eines oder mehrerer Parameter von einer Gleichung ist kein Problem.
Genau genommen liegt mein Problem darin, das lsqnonlin oder fminunc nur eine Funktion, ob in Summe von Funktionen oder nur einer Funktion,
"Auswerten" können. Ich habe ja 4 Lösungsfunktionen/-vektoren und Messdaten-Vektoren!!!!
Weiss jemand, wie man das mit MAtlab/Octave umsetzen kann?
Vielen Dank
rob
Hier mal mein Code (funktionierender Octave-Code ) für die Lösung, mit nur einer DGL. Der Code besteht aus 3 M-Files.
Code:
% M-File 1 -> wird im Promt aufgerufen
a0 = 1;
fun1 = @(a0) mkq(a0);
k = fminunc(fun1,a0);
% -- Ende M-File 1
% --M-File 2--> Funktion, aus der die Lösung für jedes a der DGL berechnet wird function lsg = mkq(a) % -- Lösung der DGL -----------
t = linspace(0,300,50);
x_start = 40;
fun2 = @(x_start,t) dgl_mkq(x_start,t,a);
f_modell = lsode(fun2,x_start,t);
% -- Ende Lösung der DGL---------------
y_mess = linspace(0,300,50) ; % Vekor der Messdaten -> fiktiv for i=1:1:length(t)
F(i) = (f_modell(i)-y_mess(i))^2;
end
lsg = sum(F);
end % -- Ende M-File 2
% --M-File 3 function dy1 = dgl_mkq(x1,t,a1)
dy1 = -a1*x1;
end
wenn lsode auch Systeme von DGLen lösen kann, ist das doch kein Problem? Du hängst die vier Lösungsvektoren untereinander und vergleichst mit den ebenso untereinandergehängten Daten.
Was meinst du genau mit "Vektoren untereinander hängen"?
Meinst du, das ich die Summe aller Fehlerquadrate aller Vektoren (fminunc) oder
alle Fehlerquadrate aller Vektoren in einem Vektor (lsqnonlin) fminunc oder lsqnonlin übergebe?
So zu sagen, das die Summe der Fehlerquadrate die Summe der Fehlerquadrate jeder einzellnen Gleichung im DGl-System ist und im "ganzen" gegen Null gehen muss, um das Problem zu Lösen?
Wenn du das meinst, darüber habe ich ach schon nachgedacht, habe aber gedacht, dass das mathematisch vll nicht korrekt ist???!!
genau das meine ich, und ich sehe darin aus mathematischer Sicht kein Problem. Die andere Frage wird sein, ob und wie schnell die Löser konvergieren, aber das wirst du ausprobieren müssen.
Hätte da noch ein kleines Problem. Umd zwar variert fminunc oder lsqnonlin die Parameter a1-5 nicht. Es werden immer die gleichen Parameterwerte (Anfangswerte) verwendet.
Das kann mann sehen wenn man sich die Mess- und Modellkurven plottet. Die Modellkurven werden einfach nicht an die Messkurven angepasst!
Ich kann mir das einfach nicht erklären! Hat jemand eine Idee, woran das liegt und was mach dagegen tun kann?
Ich hab mir dazu mal Variablen angelegt, die alle von fminunc verwendeten Parameter speichern. speicher_a1-5.
% --- Zählen aller a-Werte die von fminunc gewählt werden if(isempty(a0(1))) disp('a1 wurde nicht übergeben');
else
count_a1 = count_a1+1;
if(count_a1)>0% es wird nur gespeichert, wen a übergeben und ein hochgezählt wurde
speicher_a1(count_a1) = a0(1); % aktueller a1 wird gespeichert end end
if(isempty(a0(2))) disp('a2 wurde nicht übergeben');
else
count_a2 = count_a2+1;
if(count_a2)>0
speicher_a2(count_a2) = a0(2);
end end
if(isempty(a0(3))) disp('a3 wurde nicht übergeben');
else
count_a3 = count_a3+1;
if(count_a3)>0
speicher_a3(count_a3) = a0(3);
end end
if(isempty(a0(4))) disp('a4 wurde nicht übergeben');
else
count_a4 = count_a4+1;
if(count_a4)>0
speicher_a4(count_a4) = a0(4);
end end
if(isempty(a0(5))) disp('a5 wurde nicht übergeben');
else
count_a5 = count_a5+1;
if(count_a5)>0
speicher_a5(count_a5) = a0(5);
end end
% --- Ende Zählen
%--- Lösung der DGL -----------
t = [0103060120180240300360];
fun = @(t,x_start) dgl_mkq(t,x_start,a);
[t,f_modell] = ode45(fun,t,x_start);
%--- Ende Lösung der DGL---------------
% Im Experiment wurde für Gleichung y1(t) des System keine Messwerte genommen % -> Annahme Messdaten = Modelldaten!!!
y1_mess = f_modell(:,1);
y2_mess = [291264242216199193191193194];
y3_mess = [9.5730.426.4315.96.84.63.73.132.77];
y4_mess = [1.151.531.871.832.331.672.072.171.87];
% n -> Spalte des Lösunsgvektor aus Lösungsmatrix
n = 1;
for k=1:1:n*length(f_modell(:,n)) for i=1:1:length(f_modell(:,n))
F(k) = (f_modell(i,n)-y1_mess(i))^2;
end end
n = 2;
for k = (n-1)*length(f_modell(:,n))+1 :1: n*length(f_modell(:,n)) for i=1:1:length(f_modell(:,n))
F(k) = (f_modell(i,n)-y2_mess(i))^2;
end end
n = 3;
for k = (n-1)*length(f_modell(:,n))+1:1:n*length(f_modell(:,n)) for i=1:1:length(f_modell(:,n))
F(k) = (f_modell(i,n)-y3_mess(i))^2;
end end
n = 4;
for k = (n-1)*length(f_modell(:,n))+1:1:n*length(f_modell(:,n)) for i=1:1:length(f_modell(:,n))
F(k) = (f_modell(i,n)-y4_mess(i))^2;
end end
man sollte globale Variablen möglichst vermeiden und ausnutzen, dass MATLAB eine vektorisierte Sprache ist, d.h. man kann for-Schleifen größtenteils vermeiden.
Gegenwärtig werden die F(k) in den for-Schleifen ständig überschrieben, was wohl kaum beabsichtigt ist. Das kann zumindest eine Ursache für die Probleme sein.
Hast du dir auch mal die Textausgabe von fminunc bzw. das dritte Argument (exitflag) angesehen? Ich würde vermuten, dass lsqnonlin hier geeigneter ist.
Grüße,
Harald
Einstellungen und Berechtigungen
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