Verfasst am: 28.03.2014, 12:07
Titel: Verwendung von "lsqnonlin" zur Kalibrierung
Hallo an alle,
also ich habe ein Problem mit der Verwendung von "lsqnonlin". Mit Hilfe der zahlreichen und äußerst hilfreichen Beiträge in diesem Forum sowie natürlich der Matlab-Hilfe habe ich bereits alles versucht und muss mich nun trotzdem hilfesuchend an euch wenden.
Kurz zu meinem Problem:
Ich muss ein Programm zur automatischen Kalibrierung eines Simulationsmodells der Software EPA SWMM schreiben. Ich habe den Startvektor "P0" eines Parameters, der mit "lsqnonlin" so variiert werden soll, dass am Ende der Vektor "P" herauskommt, bei dem das Modell mit der Messung übereinstimmt. Weiters habe ich noch einen Vektor mit den Messwerten des Wasserstands im Kanal. Da ich keine Formel für den Zusammenhang zwischen Parameter und Wasserstand habe, muss ich innerhalb des Programms die Simulation laufen lassen, um Ergebnisse zu erhalten.
Also habe ich die Funktion "F = parvar" geschrieben (siehe Code). Diese Funktion sollte folgendes machen:
1. das fix abgespeicherte SWMM-Modell laden
2. die geänderten Parameter aus "P" im Modell ersetzen
3. Erstellen einer neuen Input-Datei für SWMM
4. Durchführung der Simulation mit SWMM
5. Auslesen des berechneten Wasserstands aus dem Report-File
6. Berechnung des Vektors "F" mit der Abweichung zwischen dem gemessenen und dem berechneten Wasserstand
Code:
function F = parvar(P) % P ist der Vektor mit den Zahlenwerten des Parameters
% 1. Laden der Modelldaten zum schreiben der inp-Datei load('Modell.mat');
% 2. Ersetzen der entsprechenden Werte des Parameters in den Modelldaten
Q = arrayfun(@num2str, P, 'unif', 0); % Umwandlung der Zahlen in Strings
for i=1:length(P)
Modell.SUBCATCHMENTS{i,6}=Q{i,1}; % Einsetzen der String-Werte end
% 3. Schreiben eines neuen inp-Files mit den ersetzten Parametern
inpwrite('INPneu.inp',Modell);
% 4. Durchführung der Simulation und Erstellen des neuen rpt-Files dos(['swmm5 ' 'INPneu.inp' 'INPneu.rpt']);
% 5. Auslesen des neuen rpt-Files der aktuellen Berechnung
resmodell = rptread('INPneu.rpt');
k = length(resmodell.LinkHBSM14);
WSPHBSM14Modell = zeros(k,1); % Vektor mit Simulationswerten
for i = 1:k
WSPHBSM14Modell(i) = str2num(resmodell.LinkHBSM14{i,5});
end
Durch das Minimieren der Abweichungen im Vektor "F" sollte der Vektor mit den Parametern "P" verändert werden, was bei mir aber nicht passiert. Das Programm rechnet zwar eine Zeit lang, gibt mir aber am Ende wieder genau den Startvektor "P0" aus. Auch der Fehler bzw. die Abweichung bleibt immer gleich und wird nicht minimiert.
Innerhalb der Funktion "parvar" verwende ich nur die zwei Funktionen "inpwrite" zum Schreiben des Input-Files und "rptread" zum Lesen des Report-Files, die aber 100%-ig funktionieren. Da kann der Fehler also nicht liegen.
Setze mal 'Display' auf 'iter', um dir anzeigen zu lassen, was im Detail passiert.
Schließlich kannst du noch mit dem Debugger überprüfen, ob unterschiedliche P-Werte in die Funktion übergeben werden (sollte an sich) und ob unterschiedliche F-Werte zurückgegeben werden (sollte an sich auch). Falls für verschiedene P gleiche F herauskommen, Schritt für Schritt durchgehen und schauen, ab wann die Ergebnisse gleich werden.
Es wäre zum Beispiel natürlich problematisch, wenn modell.mat auch eine Variable P enthält.
Grüße,
Harald
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 28.03.2014, 14:15
Titel:
Hallo Harald,
danke für deine rasche Rückmeldung! Du bist wirklich sehr zuverlässig!
Ja, es wird genau "P0" zurückgegeben. Für verschiedene Startvektoren "P0" erhalte ich auch verschiedene Abweichungen "F". Nur bleibt dieser Vektor während der Berechnung immer gleich, da "P" anscheinden immer gleich dem Startvektor "P0" ist. Es wird also nie ein veränderter Vektor "P" übergeben, was ich nicht verstehen kann.
Wenn ich 'Display' auf 'iter' setze kommt am Ende der Berechnung folgende Meldung:
Code:
Norm of First-order
Iteration Func-count f(x)step optimality CG-iterations
01008.719610
Optimization terminated: first-order optimality less than OPTIONS.TolFun,
and no negative/zero curvature detected in trust region model.
Es wird also nie ein veränderter Vektor "P" übergeben, was ich nicht verstehen kann.
Woraus schließt du das, wenn du nicht weißt, wie man den Debugger verwenden kann?
Zitat:
Heißt das, dass gar nie iteriert wird?
Das heißt, dass nach der Initialisierung abgebrochen wird. Setze mal 'TolFun' auf einen niedrigeren Wert, z.B. 1e-12.
Zitat:
Was muss ich eingeben, um den Debugger richtig zu verwenden?
Eingeben gar nichts. Du musst lediglich in einer Zeile, beispielsweise am Anfang von parfar, vorne auf den Strich klicken. Dann kannst du schrittweise durch die Datei navigieren. Das steht aber auch in der Hilfe, und die sollte man gerade als Nicht-Spezialist auch konsultieren, wenn man etwas nicht versteht ;)
Grüße,
Harald
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 31.03.2014, 10:25
Titel:
Hallo,
ich habe nun mit Hilfe der Breakpoints alle Schritte in der Funktion einzeln untersucht. Nach meinem Ermessen macht die Funktion genau das, was sie machen soll:
Das Input-File wird im ersten Durchlauf korrekt mit den Parametern aus dem Startvektor "P0" erstellt, die Simulation wird durchgeführt, das Report-File wird erzeugt und ausgelesen, die Wasserstände miteinander verglichen und zuletzt der Vektor "F" mit den Abweichungen berechnet.
Beim darauffolgenden Durchlauf wird als Vektor "P" allerdings wieder exakt derselbe Vektor, wie der Startvektor "P0", verwendet und alles ist wieder gleich, wie beim ersten Durchlauf des Programms, d.h. auch "F" bleibt gleich.
So läuft es dann ca. 11 min lang, wobei immer wieder mit SWMM simuliert wird, und gibt am Ende als gesuchten Vektor "P" eben wieder "P0" aus.
Auch der Tipp mit ändern des Wertes von 'TolFun' hat leider nichts verändert, egal welchen Wert ich eingeben habe.
Beim darauffolgenden Durchlauf wird als Vektor "P" allerdings wieder exakt derselbe Vektor, wie der Startvektor "P0", verwendet
Hast du format long verwendet, um die volle Stellenanzahl zu sehen?
Was passiert, wenn du einen leicht anderen Startvektor verwendest?
Grüße,
Harald
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 31.03.2014, 13:04
Titel:
Hallo,
Zitat:
Hast du format long verwendet, um die volle Stellenanzahl zu sehen?
Ja, habe ich und sie sind vollkommen ident! Da der Vektor 99 Werte beinhaltet, habe ich zur einfacheren Kontolle den Startvektor "P0" mit dem Ergebnisvektor "P" subtrahiert und einen reinen Nullvektor erhalten.
Zitat:
Was passiert, wenn du einen leicht anderen Startvektor verwendest?
Es macht keinen Unterschied, welchen Startvektor "P0" ich verwende, das Ergebnis ist immer ein zu "P0" völlig identer Ergebnisvektor "P". Dazu habe ich mir eine eigene Funktion programmiert, die mir vor Aufruf von "lsqnonlin" einen zufälligen Startvektor "P0" erzeugt, dessen Werte zufällig das 0,5- bis 2-fache des richtigen, d.h. des gesuchten Parameters sind.
% Erstellen des zufällig geänderten Startvektors des Parameters (Faktor zw. 0,5 und 2) for i=1:length(Modell.SUBCATCHMENTS)
P0(i) = str2num(Modell.SUBCATCHMENTS{i,6})*(0.5+1.5*rand);
end
P0 = transpose(P0); % Startvektor des Parameters mit zufälligen Zahlenwerten save('P0.mat','P0')% Speichern zum späteren Vergleich mit P
ich weiß nicht, ob es einen Unterschied macht, aber ich würde für die Schranken lb und ub immer Vektoren verwenden.
Ansonsten bin ich überfragt. Ich bräuchte da ein reproduzierbares Beispiel, was aber wohl mit diesem SWMM schwierig wird.
Alternativ kann vielleicht der Technische Support von MathWorks weiterhelfen.
Grüße,
Harald
Gast
Beiträge: ---
Anmeldedatum: ---
Wohnort: ---
Version: ---
Verfasst am: 31.03.2014, 14:55
Titel:
Zitat:
Ich weiß nicht, ob es einen Unterschied macht, aber ich würde für die Schranken lb und ub immer Vektoren verwenden.
Hat leider auch nichts bewirkt.
Trotzdem vielen Dank an Dich, Harald, für Deine Zeit und Dein Bemühen!
Falls ich doch noch eine Lösung oder einen anderen Weg finden sollte, werde ich es Dich und alle anderen interessierten natürlich wissen lassen.
eine Möglichkeit vielleicht noch: kannst du den SWMM-Teil vorübergehend mal durch ein anderes Modell ersetzen und schauen, ob du das Problem dort auch hast?
Wäre gut, wenn man das Problem wirklich vor sich hat und reproduzieren kann.
Versuche auch mal folgendes:
Code:
P = fminsearch(@(P)norm(parvar(P)),P0,1,800,options)% oder
P = fmincon(@(P)norm(parvar(P)),P0,1,800,options)
Verfasst am: 04.04.2014, 10:42
Titel: Problem gefunden!
Hallo nochmal,
gestern habe ich gemeinsam mit einem Matlab-Spezialisten mein Programm untersucht. Dabei hat er festgestellt, dass mein Programm und damit auch die Optimierung zwar richtig läuft und die Formulierungen alle in Ordnung sind, die ganze Sache aber dennoch einen riesen Haken hat.
Das Problem liegt in der Rechengenauigkeit des verwendeten Simulationsporgramms SWMM. Ich hoffe ich kann das richtig und verständlich erklären: Da SWMM nur mit 3 Nachkommastellen gerechnet hat, hat sich innerhalb der Funktion "lsqnonlin" bei der Berechnung der Ableitungen der Parameter zum Erstellen der Jacobimatrix nichts verändert, da die Änderungen in einem wesentlich kleineren Bereich stattfinden. Da also der Fehler somit konstant ist, war die Ableitung immer 0 und "lsqnonlin" sieht das Abbruchkriterium bereits als erfüllt an. Deshalb hat es auch gar nicht erst iteriert. Die erhaltene Meldung am Ende der Berechnung
Code:
Norm of First-order
Iteration Func-count f(x)step optimality CG-iterations
01008.719610
Optimization terminated: first-order optimality less than OPTIONS.TolFun,
and no negative/zero curvature detected in trust region model.
zeigt genau das an. Es wurden genau 100 Funktionsaufrufe zum Erstellen der Jacobimatrix (99 Parameter + 1 Funktionswert = 100 Funktionsaufrufe) durchgeführt und, da die Ableitungen immer 0 waren, die Optimierung abgebrochen.
Zur Lösung dieses Problems hätten wir daraufhin eine Benutzerdefinierte Jacobimatrix in die Funktion eingebaut, sodass es auch für weniger Nachkommastellen funktioniert. Nun läuft das Programm zwar richtig, iteriert auch einige Male und gibt am Ende einen veränderten Vektor "P" zurück, leider hat sich aber herausgestellt, dass das einfach zu ungenau ist und die Abweichung nicht ausreichend genug minimiert und damit die Parameter nicht genau genug angepasst werden.
Für mich gilt es nun irgendwie herauszufinden, wie ich aus dem Ergebnis von SWMM, d.h. im Report-File mehr Nachkommastellen herausholen kann. Ansonsten muss ich mir etwas anderes einfallen lassen.
Hier wird das Problem im Zusammenhang mit DGLen geschildert, das sollte sich aber auch auf SWMM übertragen lassen.
Grüße,
Harald
P.S.: bitte den Status des Beitrags aktualisieren. Wenn etwas auf "beantwortet" steht, dann geht man nicht davon aus, dass noch weitere Fragen bestehen.
Ich habe es bei meinem Parameteridentifikationsproblem (Simulinkmodell mit nichtlinearem dynamischen System+lsqnonlin) versucht mit deinem Tipp FinDiffRelStep, jedoch hilft es bei mir leider nicht. Matlab gibt immer noch die Startwerte zurück.
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