Verfasst am: 04.02.2013, 03:04
Titel: ode45 Solver, Werte zur Laufzeit auslesen und übergeben?
Hallo Leute,
ich möchte mit Matlab einen Regler samt Regelstrecke simulieren. Dazu habe ich die Regelstrecke, welche ein doppeltes Pendel ist, bereits programmiert. Diese löse ich mit dem ode45-solver.
Nun arbeiten die ode-Solver allerdings in der Form, dass sie gleich das gesamte System über der Zeit lösen. Also ziemlich statisch. Ich bekomme nur fertige Ergebnisse geliefert, die ich nicht mehr beeinflussen kann.
Ich möchte jedoch zur Laufzeit der Solver Messwerte, hier sind es Winkel, auslesen, mit denen dann der Regler arbeiten soll.
Der Regler soll dann entsprechend nachregeln und seine Werte wieder der Strecke, die vom Solver quasi "dauergelöst" wird, übergeben.
Habe ich vielleicht ein Verständnisproblem, was die Ode-solver betrifft?
Kann man zur Laufzeit evtl. gar nichts verändern? bzw. auslesen?
Allerdings erübrigt sich dann die Frage wie ich sonst mein DGL-System lösen und zur Laufzeit Parameter aus- bzw. übergeben kann.
Verfasst am: 04.02.2013, 08:49
Titel: Re: ode45 Solver, Werte zur Laufzeit auslesen und übergeben
Hallo Febreze,
Und ein weiteres Mal gebe ich zu bedenken, dass ODE-Solver eine glatte rechte Seite benötigen, da andernfalls der für die Qualität der Ergebnisse grundlegende Schrittweiten-Kontroller scheitern muss. Dann können 3 Effekte auftreten: 1. der Integrator bricht mit einer Fehlermeldung ab (das ist noch die wissenschaftlich brauchbarste Lösung), 2. der Integrator erkennt den Sprung und reduziert die Schrittweite so lange, bis der Rundungsfehler den Diskretisierungsfehler dominiert (das Ergebnis ist dann von fragwürdiger Genauigkeit), oder 3. der Integrator erkennt den Sprung nicht und "bügelt" einfach darüber hinweg (dann ist das Ergebnis sehr fragwürdig).
Jede Art von unstetiger Regelung sollte also vermieden werden oder man weicht auf einen Fixed-Step-Solver aus.
Das betrifft übrigens sehr viele Runge-Kutta-Integratoren, da nur die wenigsten eine aufwändige und niemals 100%ig funktionierende Unstetigkeitserkennung besitzen.
Hallo Jan,
vielen Dank für deine Antwort.
Auf eine allzu hohe Genauigkeit kommt es mir nicht unbedingt an. Es soll lediglich erst einmal geprüft werden, ob der Regler auch mit der Strecke klar kommt. Später soll die Strecke dann Hardware mäßig aufgebaut werden.
Vielleicht nochmal zum Verständnis meines Problems:
Bei der Strecke handelt es sich um ein doppeltes Pendel. Der Regler soll dieses so ausregeln, dass beide Pendelstangen senkrecht nach oben zeigen und so balanciert werden. Dies geschieht über das Messen der Winkel der Pendelstangen. Diese Hardware steht aber noch nicht zur Verfügung.
Daher habe ich die Strecke mit dem zugrunde liegenden DGL-System in Matlab programmiert und mit dem ode45-Solver gelöst. Es ist auch alles prima. Jedoch verstehe ich nicht, wie ich an einzelne diskrete Ergebniswerte, also die Winkel der Pendelstangen, während der Solver arbeitet herankommen kann.
Die ode45 Funktion löst ja sofort alles über den Zeitraum, ohne das man Werte manipulieren kann.
Da ich regelungstechnisch aber zur Laufzeit eingreifen muss, stehe ich gerade total auf dem Schlauch.
Nach deiner Aussage ist es also dem Lösungsverfahren der Solver geschuldet und somit gar nicht möglich Werte auszugeben bzw zu manipulieren während der Solver läuft. Das ist ja gerade wirklich schockierend!
Ein Doppelpendel ist instabil und chaotisch. Es mag sein, dass es Dir zunächst nicht auf eine hohe Genauigkeit ankommt. Nur hat dann die Simulation sehr schnell überhaupt keinen Bezug zum realen System mehr.
Zitat:
Jedoch verstehe ich nicht, wie ich an einzelne diskrete Ergebniswerte, also die Winkel der Pendelstangen, während der Solver arbeitet herankommen kann.
Das verstehe ich nicht. Die Winkel stehen doch in der zu integrierenden Funktion als Inputs zur Verfügung.
Zitat:
Die ode45 Funktion löst ja sofort alles über den Zeitraum, ohne das man Werte manipulieren kann.
Man kann innerhalb der zu integrierenden Funktion anstellen was immer man möchte. Allerdings muss es wegen der Spezifikationen des ODE45-Solvers glatt sein. Eine stetige Abhängigkeit der Kraft von einem Winkel wäre also kein Problem. Und wenn es unstetige Änderungen gibt, kann man die mit einer Event-Funktion abfangen und den Integrator neu starten.
Der ODE45 Solver wurde dazu konzipiert, glatte Funktionen effizient zu integrieren. Dazu wird für jeden (Zeit-)Schritt die Funktion mehrfach an Zwischenstellen ausgewertet.
Eventuell solltest Du also nicht schockiert sein, sondern einfach ein üassenderes Tool wählen, das zu Deinem Problem passt. Wie wäre ein Fixed-Step-Solver, vielleicht sogar ein einfaches Euler-Verfahren?
Hallo Jan, vielen Dank für Dein bemühtes Antworten.
Um es zuzugeben, ich kenne mich mit den Solvern aus Matlab nicht aus. Ich wüsste jetzt auch gar nicht, welcher Solver oder Lösungsverfahren mein Problem beheben würde. (Im Studium werden ja immer nur möglichst ideale und einfache Dinge behandelt und wenn man dann selber vor einem realen Problem steht, ist alles immer ganz anders... Aber ich möchte mich ja hier gar nicht beschweren...)
Was ich mir von irgend einem Solver oder Lösungsverfahren wünsche, ist "nur", dass dieser mir zum Zeitpunkt t1 die Parameter x(t1), phi1(t1), phi2(t1) gelöst in einen Vektor schreibt.
Dann kann ich drauf zugreifen und diese mit einem Regelalgorithmus auswerten.
Der Solver löst nun zum Zeitpunkt t2 => x(t2), phi1(t2), phi2(t2)
Jetzt lese ich die Werte wieder aus und gebe sie dem Regler... usw.
Also ganz einfach und sukzessiv gedacht.
Somit ist auch nicht unbedingt eine zeitliche Begrenzung gegeben.
Gibt es solch einen Solver oder Lösungsverfahren? Du hattest ja Euler angesprochen. Irgendwie muss ich ja mein System nach x(t), phi1(t), phi2(t) lösen können.
Anosten würde ich das System versuchen mit Simulink zu lösen. Nur da steckt ja auch wieder irgend ein Solver hinter.
Hier einmal mein Quelltext zum Doppelpendel. (Reibung in den Pendellagern und die Motoransteuerung für den Kraftimpuls wurden noch vernachlässigt. Es geht mir in erster Linie ums Prinzip)
% => an dieser Stelle liegt mein Verständnisproblem. % Die Option 'OutputFcn' verdeutlicht dies sehr gut. Ode45 löst alles und dann kann man erst auf alle Werte zugreifen.
% Wertübergabe der Lösungen
x_pos=x(:,1);
phi1=x(:,2);
phi2=x(:,3);
% Ausgabe figure(1) plot(t , ingrad(phi1)) , grid on , title('phi_1(t) obere Pendelstange') xlabel('t in [s]'), ylabel('\phi(t) in [°]')
figure(2) plot(t , ingrad(phi2)) , grid on , title('phi_2(t) untere Pendelstange') xlabel('t in [s]'), ylabel('\phi(t) in [°]')
Für die Verwendung komplizierter numerischer Verfahren benötigt man eine spezielle Ausbildung in Numerik. Hier im Forum findet man viele Beispiele, bei denen bereits die Addition von Gleitkomma-Zahlen zu Verwirrung führt, da bei begrenzter Genauigkeit das Kommutativ-Gesetz nicht mehr gilt. Von Integratoren oder gar Optimierern/Parameterschätzern geht dann ein weit höheres Verwirrungspotential aus.
Deswegen sollte man bei wissenschaftlichen Arbeiten mit numerischen Methoden unbedingt einen ausgebildeten Numeriker befragen. Und als solcher weise ich nun nochmal darauf hin, dass dies kein sinnvoller Input für ODE45 ist:
Statt den Integranden unstetig zu machen integriert man zunächst bis dt, und verwendet die Endposition als Startwert für das zweite Intervall ab dt, in dem nun ein neuer Kraft-Wert vorliegt. Bei einem Euler-Verfahren würde man auch keinen Schritt über den Zeitpunkt dt hinweg festlegen und erwarten, dass das Ergebnis realistisch ist.
An welcher Stelle sollen denn nun die Kräfte der Reglers wirken? Mit welcher Frequenz ändert der Regler den seinen Output?
Statt den Integranden unstetig zu machen integriert man zunächst bis dt, und verwendet die Endposition als Startwert für das zweite Intervall ab dt, in dem nun ein neuer Kraft-Wert vorliegt.
Der Satz war Gold wert! Jetzt habe ich verstanden wie man zumindest mit den Solvern arbeiten kann. Das Programm funktioniert jetzt auch so, wie ich es mir vorstelle.
Wegen des von Dir bemängelten F0 Wertes bin ich mir trotzdem noch nicht ganz sicher, wie ich damit umgehen soll. Ich übergebe jetzt den Wert F0=0 dem Integranten. Da F0 ja bei mir direkt im Integranden steht muss ich den Wert null ja zwangsweise irgendwie übergeben, da ich sonst eine konstante Kraft hätte und keinen kurzen Impuls über dt.
Später wird ja der F0-Wert ja auch die Stellgröße sein.
Ich habe noch einmal den Code angehängt vllt. kannst Du ja nochmal drüberschauen und deine Kritik dazu äußern. Ich möchte ja etwas lernen!
Ja, genau. Nun wird nicht mehr über eine Unstetigkeit hinweg integriert.
Nun noch weniger wichtige Tips:
1. Das Anhängen des Parameters als zusätzliche Inputs ist zur Rückwärtskompatibilität immer noch möglich, aber eine anonyme Funktion ist effizienter und sauberer. Also statt:
2. Ich glaube nicht, dass es in diesem Programm zu merklichen Unterschieden kommt, aber generell kann Matlab's JIT Accelerator nur Code optimieren, wenn jeweils ein Befehl pro Zeile verwendet wird. Du schreibst nur die Zuweisungen der Variablen zusammen in eine Zeile und hier gibt es kaum etwas zu beschleunigen. Falls dies innerhalb einer Schleife verwendet wird, ist der Unterschied aber messbar bis merkbar.
3. In "dxdt = x(length(x(:,1)),:);" wird der Vektor "x(:,1)" explizit erzeugt, was Speicherplatz und Rechenzeit erfordert. Eigentlich möchtest Du aber nur die Größe der ersten Dimension wissen:
Die genaue Funktion des JIT Accelerators ist nicht dokumentiert, und zwar mit Absicht. Denn TMW möchte nicht, dass Programmierer sich nach den Features der JIT richten, sondern umgekehrt die JIT an die Programme anpassen.
Als die JIT noch neu war, also z.B. in Matlab 6.5 hat der Profiler angezeigt, welche Zeilen aus welchem Grund nicht beschleunigt werden konnten. Seit dem hat sich die Beschleunigung zwar weiterentwickelt, aber die Grundzüge sind geblieben. Ich habe einige Test-Funktionen geschrieben, die dies in neuen Releases nochmal überprüfen, aber das bleibt immer etwas vage: Es bringt nichts eine 3-Zeilen-Schleife ausgiebig zu testen, wenn die JIT gerade in größeren Schleifen durch Umsortierung von Zeilen Zeit sparen kann.
Gruß, Jan
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