Dynamische Systeme in Python modellieren
Florian Lindner
ich bin momentan daran eine Anwendung zu entwickeln, die ein dynamisches
System abbilden soll. Wobei die Bezeichnung "dynamisch" erstmal nicht
wirklich korrekt ist. (da ich erstmal vor allem an stationï¿œren
Endzustï¿œnden interessiert bin)
Das System besteht aus Komponenten (Blï¿œcken) mit definierten Ein- und
Ausgï¿œngen, die mit anderen Komponenten (Serienkopplung) oder auch mit
sich selber (Rï¿œckkopplung) verbunden sind. Das Verhalten der Komponenten
ist nur abhￜngig von den Eingangsgrￜￜen, also keine Funktion der Zeit.
Ein Beispiel:
Eine Gasturbine besteht aus Verdichter, Brennkammer und Turbine.
Verdichter und Turbine sind auf einer Welle, die Turbine liefert die
Leistung um den Verdichter anzutreiben. Der Verdichter liefert (abhï¿œngig
von seiner Leistung) komprimierte Luft, die die Turbine antreibt, die
wiederum abhï¿œngig davon unterschiedliche Leistungen zur Verfï¿œgung stellt.
Verdichter zu Turbine ist also eine Kopplung ï¿œber die komprimierte Luft
und Turbine zu Verdichter eine Kopplung ï¿œber die Welle.
Irgendwo wird sich nun ein Gleichgewichtszzustand einstellen, das System
ist stationï¿œr. Leider habe ich den Zusammenhang nicht analytisch
gegeben, muss ihn also iterativ finden.
Wie ist sowas programmiere bereitet mir noch etwas Kopfzerbrechen. Ich
suche Tipps oder Hinweise auf Artikel (deutsch oder englisch) die auf so
Probleme etwas eingehen. Am Ende soll es natï¿œrlich in Python
implementiert werden, aber ich denke, das Prinzip ist fï¿œr alle Sprachen
ï¿œhnlich.
Danke,
Florian
Thomas Rachel
> Das System besteht aus Komponenten (Blöcken) mit definierten Ein- und
> Ausgängen, die mit anderen Komponenten (Serienkopplung) oder auch mit
> sich selber (Rückkopplung) verbunden sind. Das Verhalten der Komponenten
> ist nur abhängig von den Eingangsgrößen, also keine Funktion der Zeit.
Das klingt für mich wie eine Art LTI(linear time-invariant)-System, kann
das sein?
> Irgendwo wird sich nun ein Gleichgewichtszzustand einstellen, das System
> ist stationär. Leider habe ich den Zusammenhang nicht analytisch
> gegeben, muss ihn also iterativ finden.
ach so. D.h., falls Du irgendwo differentelle Zusammenhänge hast (wie
etwa Zeit/Geschwindigkeit etc., ist das nicht als DGL gegeben, sondern
als "direkter" Zusammenhang, abhängig vom alten Zustand?
Literaturhinweise hab ich da keine, aber ich würde da ungefähr so vorgehen:
* Klassen für die Subsysteme (Regler, Strecke, ...) mit einheitliches
Interface definieren (am einfachsten wohl einfach __call__
implementieren)
* Zustandsklasse definieren (entweder was NumPy-mäßiges oder eine eigene
Klasse, die Operatoren wie +,-,* definiert)
Dann kannst Du einen einfachen Regelkreis etwa so aufbauen:
setv=10 # Sollwert - evtl. auch als Funktion, falls er variabel sein soll
out=0 # oder (0,) oder (0,0) oder wie auch immer
R1=System1(...) # Regler
R2=System2(...) # Strecke
while True: # ggf. Abbruchbedingng
out=R2(R1(setv-out)))
print out # oder sonstwas mit machen
> Am Ende soll es natürlich in Python
> implementiert werden, aber ich denke, das Prinzip ist für alle Sprachen
> ähnlich.
Ja, wobei ich es mir in C etwas schwieriger vorstelle als in Python...
Thomas
Florian Lindner
> Florian Lindner schrieb:
>
>> Das System besteht aus Komponenten (Blöcken) mit definierten Ein- und
>> Ausgängen, die mit anderen Komponenten (Serienkopplung) oder auch mit
>> sich selber (Rückkopplung) verbunden sind. Das Verhalten der Komponenten
>> ist nur abhängig von den Eingangsgrößen, also keine Funktion der Zeit.
>
> Das klingt für mich wie eine Art LTI(linear time-invariant)-System, kann
> das sein?
Ich denke eher nicht. Die Beziehung zwischen Ein- und Ausgang eines
Blockes ist beliebig, also auch nichtlinear, evtl. sogar unstetig.
>> Irgendwo wird sich nun ein Gleichgewichtszzustand einstellen, das System
>> ist stationär. Leider habe ich den Zusammenhang nicht analytisch
>> gegeben, muss ihn also iterativ finden.
>
> ach so. D.h., falls Du irgendwo differentelle Zusammenhänge hast (wie
> etwa Zeit/Geschwindigkeit etc., ist das nicht als DGL gegeben, sondern
> als "direkter" Zusammenhang, abhängig vom alten Zustand?
Als DGL ist leider nichts gegeben. Man könnte wohl auch keine DGL
ableiten, da die Blöcke als Black-Boxes betrachtet werden. Prinzipiell
kann sich hinter einen Block auch ein impliziter Zusammenhang verbergen,
z.B. eine numerische Simulation einer Teilkomponente.
> Literaturhinweise hab ich da keine, aber ich würde da ungefähr so vorgehen:
>
> * Klassen für die Subsysteme (Regler, Strecke, ...) mit einheitliches
> Interface definieren (am einfachsten wohl einfach __call__
> implementieren)
Das ist schon vorhanden.
> * Zustandsklasse definieren (entweder was NumPy-mäßiges oder eine eigene
> Klasse, die Operatoren wie +,-,* definiert)
Es gibt leider nicht soetwas wie einen Zustandsvektor, die inneren
Zustände verbleiben in den jeweiligen Blöcken.
> Dann kannst Du einen einfachen Regelkreis etwa so aufbauen:
Erstmal geht es mir ja nur die Strecke und dabei darum stationäre
Gleichgewichtszustände zu entdecken.
Quasi die Blöcke, beliebig verbunden pendeln sich ausgehend von einem
Startzustand und gewissen (unveränderlichen) Umgebungsbedigungen auf
einen stabilen Zustand ein.
Ich habe PyDSTool entdeckt, mal schauen, ob das das richtige ist.
http://www.cam.cornell.edu/~rclewley/cgi-bin/moin.cgi/
Werde da mal anfragen.
> setv=10 # Sollwert - evtl. auch als Funktion, falls er variabel sein soll
> out=0 # oder (0,) oder (0,0) oder wie auch immer
> R1=System1(...) # Regler
> R2=System2(...) # Strecke
>
> while True: # ggf. Abbruchbedingng
> out=R2(R1(setv-out)))
> print out # oder sonstwas mit machen
>
>
>> Am Ende soll es natürlich in Python
>> implementiert werden, aber ich denke, das Prinzip ist für alle Sprachen
>> ähnlich.
>
>
> Ja, wobei ich es mir in C etwas schwieriger vorstelle als in Python...
Hehen, mir Sicherheit!
Grüße,
Florian
Thomas Rachel
> Thomas Rachel schrieb:
>> Florian Lindner schrieb:
>>
>>> Das System besteht aus Komponenten (Blöcken) mit definierten Ein- und
>>> Ausgängen, die mit anderen Komponenten (Serienkopplung) oder auch mit
>>> sich selber (Rückkopplung) verbunden sind. Das Verhalten der Komponenten
>>> ist nur abhängig von den Eingangsgrößen, also keine Funktion der Zeit.
>>
>> Das klingt für mich wie eine Art LTI(linear time-invariant)-System, kann
>> das sein?
>
> Ich denke eher nicht. Die Beziehung zwischen Ein- und Ausgang eines
> Blockes ist beliebig, also auch nichtlinear, evtl. sogar unstetig.
Ok. D.h. ich hab Eingang, internen Zustand und Ausgang und muß zusehen,
daß ich daraus was mache...
> Als DGL ist leider nichts gegeben. Man könnte wohl auch keine DGL
> ableiten, da die Blöcke als Black-Boxes betrachtet werden. Prinzipiell
> kann sich hinter einen Block auch ein impliziter Zusammenhang verbergen,
> z.B. eine numerische Simulation einer Teilkomponente.
Ok.
>> * Zustandsklasse definieren (entweder was NumPy-mäßiges oder eine eigene
>> Klasse, die Operatoren wie +,-,* definiert)
>
> Es gibt leider nicht soetwas wie einen Zustandsvektor, die inneren
> Zustände verbleiben in den jeweiligen Blöcken.
Die ja. Aber zwischen den Blöcken muß es ja auch was geben - evtl. auch
Summation, wie etwa die Differenzbildung bei einem gegengekoppelten
System etc.
> Erstmal geht es mir ja nur die Strecke und dabei darum stationäre
> Gleichgewichtszustände zu entdecken.
> Quasi die Blöcke, beliebig verbunden pendeln sich ausgehend von einem
> Startzustand und gewissen (unveränderlichen) Umgebungsbedigungen auf
> einen stabilen Zustand ein.
>
> Ich habe PyDSTool entdeckt, mal schauen, ob das das richtige ist.
> http://www.cam.cornell.edu/~rclewley/cgi-bin/moin.cgi/
> Werde da mal anfragen.
Ok, kenne ich noch nicht. Viel Glück!
Thomas
Florian Lindner
> Florian Lindner schrieb:
>
>> Thomas Rachel schrieb:
>>> Florian Lindner schrieb:
>>>
>>>> Das System besteht aus Komponenten (Blöcken) mit definierten Ein- und
>>>> Ausgängen, die mit anderen Komponenten (Serienkopplung) oder auch mit
>>>> sich selber (Rückkopplung) verbunden sind. Das Verhalten der Komponenten
>>>> ist nur abhängig von den Eingangsgrößen, also keine Funktion der Zeit.
>>> Das klingt für mich wie eine Art LTI(linear time-invariant)-System, kann
>>> das sein?
>> Ich denke eher nicht. Die Beziehung zwischen Ein- und Ausgang eines
>> Blockes ist beliebig, also auch nichtlinear, evtl. sogar unstetig.
>
> Ok. D.h. ich hab Eingang, internen Zustand und Ausgang und muß zusehen,
> daß ich daraus was mache...
Genau!
>> Als DGL ist leider nichts gegeben. Man könnte wohl auch keine DGL
>> ableiten, da die Blöcke als Black-Boxes betrachtet werden. Prinzipiell
>> kann sich hinter einen Block auch ein impliziter Zusammenhang verbergen,
>> z.B. eine numerische Simulation einer Teilkomponente.
>
> Ok.
>
>
>>> * Zustandsklasse definieren (entweder was NumPy-mäßiges oder eine eigene
>>> Klasse, die Operatoren wie +,-,* definiert)
>> Es gibt leider nicht soetwas wie einen Zustandsvektor, die inneren
>> Zustände verbleiben in den jeweiligen Blöcken.
>
> Die ja. Aber zwischen den Blöcken muß es ja auch was geben - evtl. auch
> Summation, wie etwa die Differenzbildung bei einem gegengekoppelten
> System etc.
Momentan plane ich diese Funktionen wiederum auch durch einzelne Blöcke
zu implementieren. Ich will ja auch kein allgemeingültiges System
erschaffen (so wie Simulink z.B.) sondern es soll passend für meine
Zwecke sein.
Ich habe da mal schnell was hingehackt:
class V(object):
T0 = 200
w = 10000
def T1(self):
return self.w + self.T0
class C(object):
T0 = 400
def T1(self):
return self.T0 + 300
class T(object):
a = 0.8
T_amb = 200
T0 = 400
def w(self):
return self.a*(self.T0-self.T_amb)
def sim():
for i in range(30):
print i
v.w = t.w()
c.T0 = v.T1()
t.T0 = c.T1()
print "T1v:", v.T1(), "T1c:", c.T1()
print "wv:", v.w, "wt:", t.w(), "delta:", abs(v.w-t.w())
print ""
v, c, t = V(), C(), T()
sim()
Es gibt also die Blöcle v, c und t. Die sind alle als serie gekoppelt
und t und v zusätlich noch über w rückgekoppelt.
Als Ergebnis kommt dann:
29
T1v: 1398.39067795 T1c: 1698.39067795
wv: 1198.39067795 wt: 1198.71254236 delta: 0.321864410214
Was auch schon ganz gut bei der analytischen Lösung liegt.
>> Erstmal geht es mir ja nur die Strecke und dabei darum stationäre
>> Gleichgewichtszustände zu entdecken.
>> Quasi die Blöcke, beliebig verbunden pendeln sich ausgehend von einem
>> Startzustand und gewissen (unveränderlichen) Umgebungsbedigungen auf
>> einen stabilen Zustand ein.
>>
>> Ich habe PyDSTool entdeckt, mal schauen, ob das das richtige ist.
>> http://www.cam.cornell.edu/~rclewley/cgi-bin/moin.cgi/
>> Werde da mal anfragen.
>
> Ok, kenne ich noch nicht. Viel Glück!
Grüße,
Florian
Florian Lindner
class DynSystem(object):
pathways = []
max_iter = 300
epsilon = 1e-4
def connect(self, input, output, *objects):
self.pathways.append( (input, output, objects) )
def run(self):
values = []
for it in range(self.max_iter):
for p in self.pathways:
objects = p[2]
for i, obj in enumerate(objects[:-1]):
val = self.get(obj, p[1]) # Value of the output
self.set(objects[i+1], p[0], val) # Set the input
values.append(val)
print "Iteration", it
if self.check_convergence(values):
break
values = []
def check_convergence(self, values):
if not hasattr(self, "old_values"):
self.old_values = values
self.old_sum_e = sys.maxint
return
sum_e = sum( map(lambda a, b: abs(a-b)**2, values,
self.old_values) )
if abs(sum_e - self.old_sum_e) < self.epsilon:
print "Convergence!"
return True
else:
self.old_values = values
self.old_sum_e = sum_e
return False
def set(self, obj, attr, value):
""" Sets the value of the attribute or function named attr from
obj. """
if attr[-2:] == "()":
getattr(obj, attr[:-2])(value)
else:
setattr(obj, attr, value)
def get(self, obj, attr):
""" Gets the value of the attribute or function named attr from
obj. """
if attr[-2:] == "()":
return getattr(obj, attr[:-2])()
else:
return getattr(obj, attr)
Benutzt wird es so:
v, c, t = V(), C(), T()
d = DynSystem()
d.connect("w", "w()", t, v) # Verbindet t zu v mit Eingang w und den
Ausgang w()
d.connect("T0", "T1()", v, c, t)
d.run()
die beiden Methoden get und set sind dafür da, zu abstrahieren, ob ich
ein Attribut setze/auslese oder eine Methode aufrufe.
Der ganze Code ist auch unter http://xgm.de/test.py verfügbar.
Freue mich über Kommentare und Verbesserungsvorschläge. So ganz
ausgereift ist er mit Sicherheit noch nicht.
Grüße,
Florian
Andreas Bruhn
Am Sat, 19 Sep 2009 13:07:55 +0200 schrieb Florian Lindner:
> Erstmal geht es mir ja nur die Strecke und dabei darum station�re
> Gleichgewichtszust�nde zu entdecken.
> Quasi die Bl�cke, beliebig verbunden pendeln sich ausgehend von einem
> Startzustand und gewissen (unver�nderlichen) Umgebungsbedigungen auf
> einen stabilen Zustand ein.
auf der Suche nach einem bestimmten Link bzw. Code hatte ich vorhin
auf meiner Platte etwas anderes gefunden. N�mlich 'SimPy'. Auszug aus
dem Readme:
:SimPy version: 1.8
:SimPy Web-site: http://simpy.sourceforge.net/
Habe noch nie mit gearbeitet, keine Ahnung ob die Version noch aktuell
ist. Aber vielleicht hilft es...
Der Link den ich eigentlich suchte, war
http://code.activestate.com/recipes/355045/
K�nnte vielleicht auch helfen. Hatte ich letztes Jahr mit der
Suchwort-Kombination 'python Spreadsheet' gefunden.
Gr��e,
Andreas