DE10205916B4

DE10205916B4 - Verfahren zur Bestimmung von Hopf-Bifurkationspunkten einer periodischen Zustandsbeschreibung eines technischen Systems

Info

Publication number: DE10205916B4
Application number: DE10205916A
Authority: DE
Inventors: Rolf Neubert
Original assignee: Infineon Technologies AG
Current assignee: Infineon Technologies AG
Priority date: 2002-02-12
Filing date: 2002-02-12
Publication date: 2004-09-16
Anticipated expiration: 2022-02-13
Also published as: US20030154225A1; DE10205916A1

Abstract

Verfahren zur Bestimmung von Hopf-Bifurkationspunkten (6) einer periodischen Zustandsbeschreibung eines Schwingungen aufweisenden technischen Systems (2), das Gleichungen der folgenden Form genügt: f(x'(t),x(r),λ) = 0dadurch gekennzeichnet, daß die folgenden Schritte ausgeführt werden:
a) Bestimmen eines Werts für den Parameter λ und Berechnen einer stationären Lösung (x,λ) durch Lösen von F(x,λ) = 0,b) Berechnen der Jacobi-Matrizen C, G,
c) Aufstellen der Funktion g ~(z) = cT(G + zC)–1b,wobei b, c orthonormierte Zufallsvektoren sind,
d) Abspalten des konstanten Terms g_∞ g ~(z) = g∞ + g(z),e) Berechnen der Windungszahl WZ der Funktion g, ausgewertet entlang der imaginären Achse:

f) Bestimmen des Monotonieverhaltens des Realteils

g(ik)) der Funktion g(ik),
g) Feststellen, ob ein Hopf-Bifurkationspunkt (6) gefunden wurde, unter Verwendung der Informationen aus den Windungszahlen WZ von g(ik) sowie unter Verwendung der Informationen aus dem Monotonieverhalten von

g(ik)): Falls Hopf-Bifurkationspunkt (6) gefunden wurde: Bestimmen einer Näherung für die Frequenz der Oszillation und gegebenenfalls...

Description

Die Erfindung betrifft ein Verfahren zur Bestimmung von Hopf-Bifurkationspunkten im Rahmen der Bestimmung einer periodischen Zustandsbeschreibung eines technischen Systems sowie dessen Verwendung. Eine Vielzahl von technischen Systemen insbesondere elektrische Schaltungen unterliegen in unterschiedlichen Zuständen Schwingungen. Solche technischen Systeme sind durch differential-algebraische Gleichungssysteme beschreibbar. Für diese technischen Systeme gilt es, sowohl stationäre als auch periodische Lösungen zu ermitteln, um eine periodische Zustandsbeschreibung des technischen Systems erstellen zu können.

Es ist ein sogenanntes QZ-Verfahren bekannt, das eine Berechnung sämtlicher Eigenwerte des zugrundeliegenden Eigenwertproblems berechnet. Dieses Verfahren ist von kubischer Komplexität in der Dimension des zugrundeliegenden technischen Systems und setzt einen sehr großen Speicherraum sowie eine sehr schnelle Rechenzeit voraus. Daher ist dieses Verfahren für große technische Systeme ungeeignet.

Weiterhin sind Verfahren bekannt, die auf iterativen Eigenwertlösern beruhen. Diese konnten bisher nicht erfolgreich zum Einsatz gebracht werden, weil sie ungenaue Lösungen liefern.

Aus Rizzoli, V.; Neri, A.: Automatic Detection of Hopf Bifurcations on the Solution Path of a Parametrized Nonlinear Circuit, In: IEEE Microwave and Guided Wave Letters. ISSN 1051-8207. 1993, Vol. 3, No. 7, S. 219-221 ist ein numerisches Verfahren für eine automatische Bestimmung von Hopf-Bifurkationen auf einem Lösungspfad eines kontinuierlich parametrisierten nicht-linearen Schaltkreises bekannt, der in einem periodisch stabilen Zustand betrieben wird. Dabei handelt es sich um eine sogenannte Torus-Verzweigung bzw. um eine sogenannte "secondary Hopf Bifurcation", also um eine Verzweigung in eine quasiperiodische Lösung von einem Zweig periodischer Lösungen.

Es ist Aufgabe der Erfindung, ein Verfahren bereitzustellen, mit dem ein Kernstück des Verfahrens zur Bestimmung periodischer Lösungen eines technischen Systems, nämlich die Bestim mung von Hopfschen Verzweigungspunkten hocheffizient und zuverlässig durchgeführt werden kann.

Diese Aufgabe wird durch den Gegenstand des unabhängigen Anspruchs gelöst. Vorteilhafte Weiterbildungen ergeben sich aus den abhängigen Ansprüchen.

Gemäß einer Voraussetzung der Erfindung ist ein Schwingungen aufweisendes technisches System durch ein System parameterabhängiger differential-algebraischer Gleichungen der Form: f(x'(t),x(t),λ) = 0gegeben. Des weiteren ist die Funktion: F(x,λ) = f(0,x,λ) = 0gegeben, welche die stationären Lösungen des Systems beschreibt.

Zunächst wird ein Bifurkationsdiagramm der Lösungen zu F(x,λ) = 0 berechnet.

Die Bestimmung eines Hopf-Bifurkationspunkts stellt einen essentiellen Teil bei der Berechnung periodischer Lösungen eines Schwingungen aufweisenden technischen Systems dar. An Hopf-Bifurkationspunkten verliert eine stationäre Lösung bzw. ein DC-Arbeitspunkt sowie ein technisches System an Stabilität, und eine stabile Oszillation setzt ein. Dabei entspringt ein Zweig periodischer Lösungen aus einem Zweig der stationären Lösung.

Numerisch werden Hopf-Bifurkationspunkte dadurch erkannt, daß ein Vorzeichenwechsel im Realteil für ein konjugiert komplexes Eigenwertpaar nachgewiesen wird. Dieses Eigenwertpaar ist aus dem Eigenwertproblem nach dem Lemma von Hopf herleitbar, das die Jacobi-Matrizen in einem stationären Punkt betrachtet. Der Wert des Imaginärteils der Eigenwerte entspricht der Frequenz der Oszillation in einer Umgebung des Hopf-Bifurkationspunkts. Aus der Eigenwertanalyse lassen sich in der Umgebung des Hopf-Bifurkationspunkts ein periodisch stationärer Zustand und durch Variation der verwendeten Parameter die Lösung im gewünschten Arbeitsbereich gewinnen.

Die Bestimmung des Bifurkationsdiagramms sowie der Hopf-Bifurkationspunkte ist aus den Dokumenten [1],[2] oder [3] bekannt.

Der nächste Schritt des erfindungsgemäßen Verfahrens sieht das Wählen eines bestimmten Wertes für den Fortsetzungsparameter λ vor. Im Nachfolgeschritt werden die Jacobi-Matrizen C, G berechnet. Diese Berechnung ist aus Dokument [4] bekannt.

Anschließend wird die Funktion: g ~(z) = cT (G + zC)–1baufgestellt. Dabei sind b und c zufällig gewählte orthonormierte Vektoren, die jedoch nicht im links- oder rechtseitigen Kern von C liegen dürfen. g ~(z) ist eine gebrochen rationale Funktion mit reellen Koeffizienten.

Von dieser Funktion g ~(z) wird der konstante Term g_∞ abgespalten. Dies erfolgt mittels einer Polynomdivision. Somit ergibt sich: g ~(z) = g∞ + g(z)

Diese Abspaltung ist dem Fachmann bekannt. Nachfolgend wird g(z)anstatt g ~(z) betrachtet. g(z) hat die gleichen Polstellen wie g ~(z) und ebenfalls reelle Koeffizienten, verschwindet allerdings im Unendlichen. Durch parallele Auswertung von g(z) mit verschiedenen, festen Zufallsvektoren b, c können diejenigen Informationen über die Pole herausgefiltert werden, die Aufschluß über die Eigenwerte geben. Eigenwerte sind ein Maß für die Stabilität und sind aus den Nullstellen des Nenners berechenbar. Es erfolgt eine Auswertung von g entlang der imaginären Achse, wobei aus Stetigkeitsgründen angenommen wird, daß kein Eigenvektor einen Realteil mit dem Wert 0 aufweist, d.h. man berechnet g(ik) fortlaufend für alle Vektoren b, c parallel, beginnend k₀ = 0 für wachsende Werte von k ∊ R solange, bis ⎪g(ik)⎪ hinreichend klein wird. i bezeichnet hier die imaginäre Einheit. Dabei wird die Windungszahl WZ sowie das Monotonieverhalten der Funktion g(ik). betrachtet.

Die Windungsfunktion W der Funktion g(ik) ergibt sich wie folgt: W(g(ik);ko,k1) = arg (g(ik1)) – arg(g(ik0))

Bei der Funktion: arg handelt es sich um die Winkelfunktion. Im Falle einer Konvergenz ergibt sich die Windungszahl WZ, d.h. der Grenzwert der Windungsfunktion W wie folgt:

Aus der Mathematik läßt sich mittels des Gaußschen Integralsatzes ein Zusammenhang zwischen der Windungszahl WZ und der Anzahl der Nullstellen und der Polstellen herleiten. Es ergibt sich:

Dabei bezeichnen N_l, N_r, P_l und P_r die Anzahl der Null- und Polstellen in der rechten bzw. linken Halbebene ihren Vielfachheiten entsprechend gezählt.

Dementsprechend erfolgt eine Berechnung der Windungszahl WZ der Funktion g(ik) bei festen Zufallsvektoren b, c zu verschiedenen Fortsetzungsparametern λ. Als Ergebnis erhält man eine Übersicht über das Verhalten der Null- und Polstellen von g(ik) während einer Fortsetzung.

Überquert ein konjugiert komplexes Polstellenpaar die imaginäre Achse, ändert sich die Windungszahl WZ konsequenterweise um ±2π. Dementsprechend können Hopf-Bifurkationspunkte durch eine Analyse der Windungszahl WZ im Laufe einer Fortsetzung, d.h. für wachsende Fortsetzungsparameter λ identifiziert werden.

In einem weiteren Schritt des erfindungsgemäßen Verfahrens erfolgt die Bestimmung des Monotonieverhaltens des Realteils

g(ik)) der Funktion g(ik).

Hierbei wird der Umstand genutzt, daß das Überqueren der imaginären Achse durch ein konjugiert komplexes Eigenwertpaar verursacht, daß sich ein lokales Minimum des Realteils

g(ik)) der Funktion g(ik) in ein lokales Maximum in einer Umgebung des Hopf-Bifurkationspunkts verändert. Dementsprechend gibt eine Untersuchung des Realteils

g(ik)) in Abhängigkeit von variierenden Fortsetzungsparametern λ Aufschluß darüber, ob ein Hopf-Bifurkationspunkt verpaßt wurde. Ein Hopf-Bifurkationspunkt liegt insbesondere dann vor, wenn die Steigung bzw. die erste Ableitung des Realteils

g(ik)) der Funktion g(ik) für aufeinanderfolgende Parameterwerte von λ einen Vorzeichenwechsel aufweist.

Während der Ausführung des erfindungsgemäßen Verfahrens werden die Windungszahl WZ sowie die Monotoniewechsel des Monotonieverhaltens des Realteils

g(ik)) für variierende Fortsetzungsparameter λ protokolliert. Unter Verwendung dieser Informationen wird in dem nun folgenden Schritt festgestellt, ob ein Hopf-Bifurkationspunkt gefunden wurde. Hierbei wird auf die Informationen aus den Windungszahlen WZ und deren Zusammenhang zwischen den Nullstellen N und den Polstellen P von g(ik) sowie auf die Informationen über das Monotonieverhalten von

g(ik)) zugegriffen.

Falls ein Hopf-Bifurkationspunkt gefunden wurde, wird die Frequenz der Oszillation in der Umgebung des Hopfschen Verzweigungspunktes bestimmt. Das erfindungsgemäße Verfahren ist mit der Bestimmung des Hopf-Bifurkationspunkts an dieser Stelle beendet.

Falls kein Hopf-Bifurkationspunkt gefunden wurde, erfolgt eine Überprüfung, ob der Parameter λ einen bestimmten durch einen Benutzer vorgebbaren Wert λ_max erreicht hat. Ist dies der Fall, so ist das erfindungsgemäße Verfahren ebenfalls beendet.

Falls kein Hopf-Bifurkationspunkt gefunden wurde und falls λ den Wert λ_max noch nicht erreicht hat, wird das erfindungsgemäße Verfahren ab dem Schritt des Bestimmens eines Werts für den Parameter λ wiederholt.

Ein Grundgedanke der Erfindung besteht darin, die Funktion g(ik) auszuwerten, anstelle eine Berechnung sämtlicher Eigenwerte vorzunehmen.

Ein weiterer Grundgedanke der Erfindung besteht darin, das Verfahren zur Bestimmung von Hopf-Bifurkationspunkten durch Windungszahlanalyse sowie das Verfahren zur Bestimmung von Hopf-Bifurkationspunkten durch Monotoniebetrachtungen zu einem gemeinsamen Verfahren zusammenzuführen, und mit Hilfe der Synthese beider Verfahren eine hocheffiziente und robuste Detektion von Hopf-Bifurkationspunkten zur Verfügung zu stellen.

Die erfindungsgemäße Auswertung der Funktion g ~(z) = c^T (G + zC)^–1b entspricht der Lösung eines linearen Gleichungssystems. Somit weist diese Auswertung für die im allgemeinen sehr dünn besetzten Jacobi-Matrizen G und C des untersuchten technischen Systems eine nur schwach superlineare Komplexität auf.

Das erfindungsgemäße Verfahren zur Bestimmung von Hopf-Bifurkationspunkten weist dementsprechend eine Komplexität von O(n^α), α≈1 auf. Bei den bekannten Verfahren beträgt die Komplexität O(n³). Dementsprechend arbeitet das erfindungsgemäße Verfahren ab n≈200 besonders effizient, bei n≈1000 liegen die Speedups schon bei einem Wert von ca. 100. Somit ist das erfindungsgemäße Verfahren auch für größere Schaltungen vorteilhaft anwendbar, bei denen die bekannten QZ-Verfahren versagen.

Gemäß einer Ausführungsform der Erfindung wird das Monotonieverhalten des Realteils

g(ik)) der Funktion g(ik) durch eine Extremwertuntersuchung von

g(ik)) für variierende Werte des Fortsetzungsparameters λ bestimmt.

Überquert zwischen zwei Fortsetzungsparametern λ₀ und λ₁ ein konjugiert komplexes Eigenwertpaar die imaginäre Achse, d.h. zwischen λ₀ und λ₁ liegt ein Hopfscher Verzweigungspunkt, so verändert sich ein lokales Minimum bzw. Maximum des Realteils

g(ik) der Funktion g(ik) am Parameterwert λ₀ in ein lokales Maximum bzw. Minimum am Parameterwert λ₁. Folglich dient auch dieses Krierium zur Charakterisierung von Hopfschen Verzweigungspunkten. Gleichbedeutend, jedoch leichter zu überprüfen ist das Kriterium, daß sich lokal in der Umgebung des Ex tremwertes das Vorzeichen der Ableitung des Realteils

g(ik)) der Funktion g(ik) verändert.

Die Untersuchung der Funktion g(ik) liefert folglich zwei Kriterien zur Detektion eines Hopf-Bifurkationspunktes: Zum einen eine Veränderung der Windungszahl um ±2π, zum anderen eine Veränderung des Extremwertverhaltens. Durch eine Kombination beider Kriterien wird ein sicheres und effizientes Auffinden von Hopf-Bifurkationspunkten bereitgestellt.

Gemäß einer weiteren Ausführungsform der Erfindung wird die Frequenz der Oszillation in der Umgebung eines Hopf-Bifurkationspunktes bestimmt. Als Näherung für die Frequenz wird die Stelle k_m gewählt, an welcher der Realteil

g(ik)) seinen steilsten Gradienten hat. Hierzu wird die Steigung der Sekanten von

g(ik)) durch die Punkte

auf ihr Betragsmaximum untersucht. Ausgehend von dieser Näherung wird der tatsächliche Eigenwert mittels einer inversen Iteration berechnet. Dieses Verfahren ist dem Fachmann bekannt. Aufbauend darauf können die aus dem Dokument [3] bekannten Methoden zur Berechnung von periodischen Lösungen angewendet werden.

Durch diese Ausführungsform der Erfindung kann ausgehend von einer Näherung eine periodische Lösung in der Umgebung eines Hopf-Bifurkationspunktes zuverlässig und genau bestimmt werden.

Eine weitere Ausführungsform der Erfindung sieht vor, daß das erfindungsgemäße Verfahren rückwirkend für Zwischenwerte durchgeführt wird, wenn keine sichere Aussage über das Vorhandensein eines Hopf-Bifurkationspunkts getroffen werden kann.

Dabei werden die Verfahrensschritte ab dem Bestimmen eines Werts für den Parameter λ wiederholt. Für den Parameter λ wird ein Wert gewählt, der zwischen dem letzten gewählten und dem vorletzten gewählten Wert des Parameters λ liegt. Der neue Parameter λ ist auch als arithmetisches Mittel aus dem letzten gewählten und aus dem vorletzten gewählten Wert des Parameters λ vorsehbar. Ausgehend von dem neu gewählten Wert für den Parameter λ wird nach beiden Seiten gerechnet, und zwar in Richtung des letzten gewählten Werts des Parameters λ sowie in Richtung des vorletzten gewählten Werts des Parameters λ. Des weiteren ist eine individuelle Anpassung der Abtastung bei der Auswertung der Funktion g(ik) möglich. Dementsprechend ist eine besonders zuverlässige und besonders feine Bestimmung eines Hopf-Bifurkationspunkts möglich.

Außerdem betrifft die Erfindung ein Verfahren, bei dem das zugrundeliegende technische System wenigstens eine elektrische Schaltung aufweist oder als elektrische Schaltung insbesondere als ein autonomer elektronischer Schaltkreis vorliegt.

Weiterhin betrifft die Erfindung die Verwendung des erfindungsgemäßen Verfahrens nach einem der vorhergehenden Ansprüche zur elektrischen Schaltungssimulation.

Die Erfindung ist auch in einem Computerprogramm zur Ausführung eines Verfahrens zur Bestimmung von Hopf-Bifurkationspunkten einer periodischen Zustandsbeschreibung eines techni schen Systems verwirklicht. Das Computerprogramm ist dabei so ausgebildet, daß nach Eingabe des technischen Systems durch ein System von differential-algebraischen Gleichungen ein Verfahren gemäß einem der erfindungsgemäßen Ansprüche verwirklicht ist. Dabei sind als Ergebnis des Verfahrens Hopf-Bifurkationspunkte einer periodischen Zustandsbeschreibung des technischen Systems ausgebbar. Durch diese Hopf-Bifurkationspunkte können sehr vorteilhaft Aussagen über das zugrundeliegende technische System getroffen werden.

Durch das erfindungsgemäß verbesserte Computerprogramm ergeben sich eine zuverlässige und vollständige Ermittlung der Hopf-Bifurkationspunkte sowie eine Laufzeitverbesserung gegenüber den bekannten Verfahren zur Bestimmung von Hopf-Bifurkationspunkten von technischen Systemen.

Die Erfindung betrifft außerdem ein Computerprogramm, das auf einem Speichermedium enthalten ist, das in einem Computerspeicher abgelegt, das in einem Direktzugriffsspeicher enthalten ist oder das auf einem elektrischen Trägersignal übertragen wird. Weiterhin betrifft die Erfindung einen Datenträger mit einem solchen Computerprogramm sowie ein Verfahren, bei dem ein solches Computerprogramm aus einem elektronischen Datennetz, wie beispielsweise aus dem Internet, auf einen an das Datennetz angeschlossenen Computer heruntergeladen wird.

Die Erfindung ist in den Figuren anhand eines Ausführungsbeispiels näher veranschaulicht.
1 zeigt eine schematische Darstellung eines Flußdiagramms 1 des Verfahrens zur Bestimmung von Hopf-Bifurkationspunkten einer periodischen Zustandsbeschreibung eines technischen Systems,
2 zeigt ein Schaltbild eines vollintegrierten spannungsgesteuerten Oszillators 2 gemäß einem Ausführungsbeispiel,
3 zeigt ein Bifurkationsdiagramm 3 des differentialalgebraischen Gleichungssystems, das den in 2 gezeigten spannungsgesteuerten Oszillator 2 beschreibt, gemäß dem Ausführungsbeispiel,
4 zeigt eine Realteilverlaufsdarstellung 8 einer Funktion g(ik) mit einem ersten Realteilverlauf 9, mit einem zweiten Realteilverlauf 10 und mit einem dritten Realteilverlauf 11 des in 2 gezeigten spannungsgesteuerten Oszillators 2 gemäß dem Ausführungsbeispiel.
1 zeigt eine schematische Darstellung eines Flußdiagramms 1 des Verfahrens zur Bestimmung von Hopf-Bifurkationspunkten einer periodischen Zustandsbeschreibung eines technischen Systems.
Das erfindungsgemäße Verfahren geht aus von einem System parameterabhängiger differential-algebraischer Gleichungen der Form f(x'(t),x(t),λ) = 0. Dabei werden zur Bestimmung von Hopf-Bifurkationspunkten die Lösungen der algebraischen Gleichung F(x,λ) = f(0,x,λ) = 0 untersucht.
Die Lösungen von F(x,λ) = 0 für einen variierenden Parameter λ untersucht man mit sogenannten Fortsetzungsverfahren. Die Durchführung von Fortsetzungsverfahren ist aus den Dokumenten [1],[2] oder [3] bekannt.
Im folgenden werden die Aktionen während eines Fortsetzungsschrittes, also für festes λ beschrieben. Nach Bestimmung des Wertes für λ wird die zugehörige statioinäre Lösung (λ,x) mittels Lösen von F(x,λ) = 0 bestimmt.
Daraufhin wird die Linearisierung um die Lösung, also die Jacobi-Matrizen C = f_x'(0,x,λ) und G = f_x(0,x,λ) berechnet . Diese Berechnung ist aus Dokument [4] bekannt.
In einem weiteren Schritt des Verfahrens wird eine Funktion g ~(z) = c^T(G + zC)^–1b aufgestellt . Dabei sind b und c orthonormierte Zufallsvektoren. g ~(z) ist eine gebrochen rationale Funktion mit reellen Koeffizienten.
Bei der Berechnung von periodischen Lösungen für große Systeme sind die Matrizen gewöhnlich nur schwach besetzt. Somit weist der Aufwand für die Auswertung von g ~(z) im allgemeinen eine nur schwach superlineare Komplexität auf. g ~(z) hat die Eigenwerte der Gleichung G(λ0)υ(λ0) + μ(λ0)C(λ0)υ(λ0) = 0als Polstellen. Dies ergibt sich aus der Cramerschen Regel, die in Dokument [5] beschrieben ist.
Der Zählergrad von g ~(z) hängt vom Rang und von der Struktur der Jacobi-Matrix C ab. Dies folgt aus der adjunkten Matrix und dem Laplaceschen Entwicklungssatz, der in Dokument [5] beschrieben ist.
Die Nullstellen hängen im Gegensatz zu den Polstellen von der Wahl der orthonormierten Zufallsvektoren b, c ab. Die Jacobi- Matrix C hat nur dann vollen Rang, wenn das System parameterabhängiger differential-algebraischer Gleichungen: f(x'(t),x(t),λ) = 0eine implizite, gewöhnliche Differentialgleichung ist. Bei vielen als Schaltkreisen vorliegenden technischen Systemen ist dies nicht der Fall. Hier ist die Jacobi-Matrix C normalerweise singulär. Deshalb ist die Annahme zulässig, daß Zählergrad und Nennergrad gleich sind.
Dementsprechend wird mittels einer Polynomdivision der konstante Term g_∞ von der Funktion g ~(z) abgespalten: g ~(z)=g∞ + g(z)
Diese Abspaltung ist dem Fachmann bekannt.
Dabei entspricht der Grenzwert g_∞ = limg ~(z) für z → ∞ dem Grenzwert von g ~(z) gegen unendlich. Im folgenden wird g(z) anstatt g ~(z) betrachtet. g(z) weist die gleichen Polstellen auf wie g ~(z) und hat ebenfalls reelle Koeffizienten. Allerdings verschwindet g(z) im Unendlichen.
Im nächsten Schritt des erfindungsgemäßen Verfahrens wird die Windungszahl WZ der Funktion g entlang der imaginären Achse ausgewertet, also g(ik) berechnet. Aus der Berechnung der Windungszahl WZ der Funktion g(ik) bei festen Zufallsvektoren b, c zu verschiedenen Fortsetzungsparametern λ erhält man eine Übersicht über das Verhalten der Nullstellen und der Polstel len von g(ik) und damit über das Verhalten der Eigenwerte während einer Fortsetzung.
Aus der Erkenntnis, ob sich das zugrundeliegende System vom stabilen in den instabilen Zustand bewegt oder von dem instabilen in den stabilen Zustand bewegt, kann man doppelte Nullstellen und Polstellen unterscheiden. Aus einer parallelen Berechnung von Windungszahlen WZ zu verschiedenen Vektoren b, c ergeben sich weitere Erkenntnisse über die Polstellen und über die Nullstellen. So können Polstellen und Nullstellen unterschieden werden und Hopf-Bifurkationspunkte identifiziert werden.
Die Windungsfunktion W ist wie folgt definiert: W(g(ik);ko + k1) = arg(g(ki1) – arg(g(ik0))
Im Falle einer Konvergenz ergibt sich die Windungszahl WZ zu:
Unter der ersten Voraussetzung, daß mit g eine rationale Funktion gegeben ist, deren Zählergrad kleiner ist als der Nennergrad und unter der zweiten Voraussetzung, daß keine Nullstelle N und keine Polstelle P von g einen Realteil mit dem Wert 0 aufweist, ergibt sich der folgende Zusammenhang zwischen der Windungszahl WZ und der Anzahl von Null- bzw. Polstellen:
Dabei bezeichnen N_l, N_r, P_l und P_r die Anzahl der Nullstellen N und der Polstellen P in der rechten bzw. der linken Halbebene, ihren Vielfachheiten entsprechend gezählt.
Des weiteren ist zu beachten, daß bei der Bestimmung der Windungszahl WZ die Änderungen immer als ganzzahlige Vielfache von π vorliegen. Das Überqueren der imaginären Achse durch ein konjugiert komplexes Eigenwertpaar oder durch ein konjugiert komplexes Nullstellenpaar wird durch eine Windungszahländerung von ±2π erkannt. Das Überqueren der imaginären Achse durch einfache Nullstellen N oder durch einfache Polstellen P wird durch eine Änderung der Windungszahl WZ von ±π angezeigt. Andere Differenzen der Windungszahlen WZ lassen sich durch die Kombination dieser drei Fälle erklären.
Das Überqueren der imaginären Achse durch ein konjugiert komplexes Eigenwertpaar oder durch ein konjugiert komplexes Nullstellenpaar weist auf das Vorhandensein eines Hopf-Bifurkationspunkts hin.
Aufgrund von zu großen Fortsetzungsschritten in λ kann es zu sich kompensierenden Überquerungen der imaginären Achse von Null- und Polstellen kommen, so daß das Windungszahl-Kriterium nicht anwendbar ist. Aus diesem Grund wird zusätzlich zur Berechnung der Windungszahlen das Monotonieverhalten des Realteils
(g(ik)) der Funktion g(ik) in der Nähe von lokalen Extrema untersucht.
Überquert ein konjugiert komplexes Eigenwertpaar zwischen zwei Fortsetzungsparametern λ₀ und λ₁ die imaginäre Achse, d.h. liegt ein Hopfscher Verzweigungspunkt zwischen λ₀ und λ₁, so verändert sich ein lokales Minimum bzw. Maximum des Realteils
g(ik)) der Funktion g(ik) am Parameterwert λ₀ in ein lokales Maximum bzw. Minimum am Parameterwert λ₁. Folglich dient auch dieses Kriterium zur Charakterisierung von Hopfschen Verzweigungspunkten. Gleichbedeutend, jedoch leichter zu überprüfen ist das Kriterium, daß sich lokal in der Umgebung des Extremwertes das Vorzeichen der Ableitung des Realteils
g(ik)) der Funktion g(ik) verändert.
Basierend auf den Kriterien der Windungszahl WZ und des Extremwertverhaltens wird nun überprüft, ob ein Hopf-Bifurkationspunkt gefunden wurde.
Falls diese Überprüfung ein negatives Ergebnis aufweist, erfolgt eine Überprüfung, ob der gewählte Parameter λ einen durch einen Benutzer vorgebbaren Wert λ_max erreicht hat.
Hat der Parameter λ den Wert λ_max noch nicht erreicht, erfolgt eine Wiederholung sämtlicher Verfahrensschritte ab dem Schritt des Wählens eines neuen Parameters λ.
Überschreitet der gewählte Parameter λ den vorgegebenen Wert λ_max, so ist das erfindungsgemäße Verfahren beendet. In diesem Fall wurde kein Hopf-Bifurkationspunkt detektiert.
Im Falle, daß die Überprüfung, ob ein Hopf-Bifurkationspunkt gewunden wurde, ein positives Ergebnis ergibt, so wird eine Näherung der Frequenz der Oszillation in der Umgebung des Hopf-Bifurkationspunktes bestimmt und gegebenenfalls eine inverse Iteration zur Bestimmung des exakten Werts für die Frequenz durchgeführt.
Kann an dieser Stelle keine sichere Aussage über das Vorhandensein eines Hopf-Bifurkationspunkts getroffen werden, wird das erfindungsgemäße Verfahren rückwirkend für Zwischenwerte des Fortsetzungsparameters λ durchgeführt.
Für die inverse Iteration wird als Näherung für die Frequenz der Oszillation die Stelle k_m gewählt, an der
g(ik)) den steilsten Gradienten hat.
Hierzu wird die Steigung der Sekanten von
g(ik)) durch die Punkte: x(km, g(ikm)), x(km–1, g(ikm–1)),also der Term
auf ein Betragsmaximum untersucht. Ausgehend von dieser Näherung wird der tatsächliche Eigenwert mittels einer inversen Iteration berechnet. Aufbauend darauf können mit den Verfahren aus dem Dokument [3] periodische Lösungen des zugrundeliegenden technischen Systems berechnet werden.
2 zeigt ein Schaltbild eines vollintegrierten spannungsgesteuerten Oszillators 2 gemäß einem Ausführungsbeispiel.
Der in 2 gezeigte spannungsgesteuerte Oszillator 2 stellt ein System dar, das periodisch wiederkehrende Zustände annimmt und wird zur Erzeugung von Schwingungen insbesondere von Hochfrequenzschwingungen im Mobilfunk eingesetzt.
Mit Einbeziehung parasitärer Elemente ist der spannungsgesteuerte Oszillator 2 durch ein System von 905 differentialalgebraischen Gleichungen der Form: f(x'(t),x(t),λ) = 0beschreibbar. Von einer Darstellung dieses differentialalgebraischen Gleichungssystems wird aus Gründen der Anschaulichkeit abgesehen.
Eine genaue Beschreibung des Aufbaus und der Funktionsweise des spannungsgesteuerten Oszillators 2 ist in Dokument [6] enthalten.
3 zeigt ein Bifurkationsdiagramm 3 des differentialalgebraischen Gleichungssystems, das den in 2 gezeigten spannungsgesteuerten Oszillator 2 beschreibt, gemäß dem Ausführungsbeispiel.
Das Bifurkationsdiagramm 3 weist stabile stationäre Lösungen 4, instabile Lösungen 5, einen Hopf-Bifurkationspunkt 6 und periodische Lösungen 7 auf.
In dem in 3 dargestellten Bifurkationsdiagramm 3 wird zu jedem Parameter λ eine eindimensionale Projektion der zugehörigen stationären bzw. periodischen Lösungen x des Gleichungssystems gezeichnet, das den spannungsgesteuerten Oszillator 2 beschreibt.
Die durchgezogen dargestellte Kurve symbolisiert die Menge der stabilen stationären Lösungen 4 des zugrundeliegenden differentia lalgebraischen Gleichungsystems des spannungsgesteuer ten Oszillators 2. Die gestrichelt dargestellte Linie stellt die instabilen stationären Lösungen 5 des zugrundeliegenden differential-algebraischen Gleichungsystems des spannungsgesteuerten Oszillators 2 dar. Die periodischen Lösungen 7 sind mittels Kreuzen dargestellt. Bei dem Punkt, an dem der Zweig der periodischen Lösungen 7 aus dem Zweig der stabilen stationären Lösungen 4 entspringt, handelt es sich um den Hopf-Bifurkationspunkt 6. Weiterhin geht an dem Hopf-Bifurkationspunkt 6 der Zweig der stabilen stationären Lösungen 4 in den Zweig der instabilen stationären Lösungen 5 über.
Bei dem in 3 dargestellten Bifurkationsdiagramm 3 wird die lokale Dynamik der Lösungen der algebraischen Gleichung: F(x,λ) = f(0,x,λ) = 0untersucht, die sich aus dem parameterabhängigen differentialalgebraischen Gleichungssystem: f(x'(r),x(i),λ)=0ergibt. Die Lösungen von F(x,λ) = 0 liegen als eindimensionale Untermannigfaltigkeiten vor und weisen oft eine komplizierte Struktur auf. Beispielsweise können zu einem Parameter λ mehrere Lösungen gehören, die Funktion F(x,λ) = 0 kann Singularitäten und Verzweigungen aufweisen.
Im allgemeinen ergibt sich aus der Berechnung der Menge der stationären Lösungen der Funktion F(x,λ) = 0 noch keine Aussage über das dynamische Verhalten des zugrundeliegenden Systems in einer Umgebung einer stationären Lösung. Dementsprechend erfolgt für jede stationäre Lösung eine Stabilitätsanalyse.
Die Menge der stationären Lösungen der Funktion F(x,λ) = 0 werden in Verbindung mit Informationen zu der Stabilität der einzelnen stationären Lösungen gewöhnlich in einem Bifurkationsdiagramm 3 dargestellt. Somit sind in dem Bifurkationsdiagramm 3 diejenigen kritischen Bereiche des Parameters λ erkennbar, bei denen ein Gewinn bzw. ein Verlust der Stabilität der jeweiligen stationären Lösung der Funktion F(x,λ) = 0 gegeben ist.
In dieser Erfindung wird besonders der Stabilitätsverlust an sogenannten Hopf-Bifurkationspunkten 6 betrachtet, an denen aus dem Zweig einer stationären Lösung ein Zweig periodischer Lösungen entspringt.
Die Ausführung der numerischen Berechnung des in 3 dargestellten Hopf-Bifurkationspunkts 6 erfolgt in diesem Ausführungsbeispiel gemäß der Vorgehensweise des erfindungsgemäßen Verfahrens.
Für das Differentialgleichungssystem, das den spannungsgesteuerten Oszillator 2 beschreibt, werden in jedem Schritt des Fortsetzungsverfahrens die nachfolgenden erfindungsgemäßen Verfahrensschritte durchgeführt, wobei von dem erstellten Bifurkationsdiagramm 3 ausgegangen wird:

– Bestimmen eines Werts für den Parameter λ,
– Bestimmen einer stationären Lösung (x,λ) durch Lösen von F(x,λ) = 0,
– Berechnen der Jacobi-Matrizen C, G,
– Aufstellen der Funktion g ~(z) = c^T(G + zC)^–1b,
– Abspalten des konstanten Terms g_∞,
– Berechnen der Windungszahl WZ der Funktion g, ausgewertet entlang der imaginären Achse
– Bestimmen des Monotonieverhaltens des Realteils
g(ik)) der Funktion g(ik),
– Auswerten der Kriterien über Windungszahl WZ und Extremwertverhalten, damit Feststellen, ob ein Hopf-Bifurkationspunkt 6 gefunden wurde,
– Falls Hopf-Bifurkationspunkt 6 gefunden wurde, Bestimmen einer Näherung für die Frequenz der Oszillation und gegebenenfalls Durchführen einer inversen Iteration,
– Falls kein Hopf-Bifurkationspunkt 6 gefunden wurde, Überprüfen ob λ ≤ λ_max,
– Falls kein Hopf-Bifurkationspunkt 6 gefunden wurde und falls λ ≤ λ_max, Wiederholen der Schritte ab dem Schritt des Bestimmens eines Werts für den Parameter λ für den nächsten Fortsetzungsschritt.

Auf eine Darstellung der konkreten Ausführung der einzelnen Verfahrensschritte für das 905 differential-algebraische Gleichungen umfassende differential-algebraische Gleichungssystems, das den spannungsgesteuerten Oszillator 2 beschreibt, wird aus Gründen der Anschaulichkeit verzichtet.
Das Ergebnis des erfindungsgemäßen Verfahrens, und zwar der Hopf-Bifurkationspunkt 6 ist in 3 graphisch dargestellt.
Das Nachvollziehen dieser Verfahrensschritte zur Bestimmung von Hopf-Bifurkationspunkten für das differential-algebraische Gleichungssystem, das den spannungsgesteuerten Oszillator 2 beschreibt und insgesamt 905 differential-algebraische Gleichungen umfaßt, ist für den Fachmann anhand der in diesem Dokument enthaltenen Informationen sowie anhand der zitierten Dokumente möglich.
Numerisch werden Hopf-Bifurkationspunkte 6 dadurch erkannt, daß ein Vorzeichenwechsel im Realteil für ein konjugiert komplexes Eigenwertpaar nachgewiesen wird, wobei bisher die Eigenwerte mit dem QZ-Verfahren berechnet wurden.
Die Durchführung des QZ-Verfahrens ist aus Dokument [8] bekannt. Dieses benötigt ungefähr 30n³ Operationen pro Eigenwertproblem, ist also ein Verfahren kubischer Komplexität. Für die Bestimmung des Hopf-Bifurkationspunkts 6 ist auch der im Dokument [7] enthaltene Satz anwendbar.
Eine effiziente Alternative mittels einer Windungszahl-Analyse ist für gewöhnliche Differentialgleichungen der Form x' = f(x,λ) in Dokument [9] beschrieben.
Durch das Verfahren zur Erkennung von Hopf-Bifurkationspunkten 6 durch Windungswinkel-Analyse werden häufig Hopf-Bifurkationspunkte 6 nicht erkannt, da sich die Sprünge der Windungszahl WZ aufheben. Dies ist besonders der Fall, wenn der Fortsetzungsparameter λ mit einer zu großen Schrittweite abgetastet wird. Dabei kommt es vor, daß in einem Schritt sowohl ein Eigenpaar von links nach rechts als auch ein Nullstellenpaar von rechts nach links wandert.
Durch die Verwendung des erfindungsgemäßen Verfahrens wird zur Berechnung einer Lösung des den spannungsgesteuerten Oszillator 2 beschreibenden differential-algebraische Gleichungssy tems eine "CPU"-Zeit bzw. eine "Central Processing Unit"-Zeit von 23.77 Sekunden benötigt. Das bekannte QZ-Verfahren benötigt zur Berechnung einer Lösung des gleichen Differentialgleichungssystems eine Zeitspanne von 1977 sec.
4 zeigt eine Realteilverlaufsdarstellung 8 einer Funktion
g(ik)) mit einem ersten Realteilverlauf 9, mit einem zweiten Realteilverlauf 10 und mit einem dritten Realteilverlauf 11 des in 2 gezeigten spannungsgesteuerten Oszillators 2 gemäß dem Ausführungsbeispiel.
Mittels der in 4 grafisch dargestellten Untersuchung des Realteils
(ik)) der Funktion g(ik) kann eine Aussage darüber getroffen werden, ob ein Hopf-Bifurkationspunkt 6 verpaßt wurde. Diese Untersuchung des Realteils von g(ik) basiert auf dem folgenden Satz:
Das Überqueren der imaginären Achse durch ein konjugiert komplexes Eigenwertpaar verursacht, daß sich ein lokales Minimum (Maximum) vom Realteil der Transferfunktion g(ik) in ein lokales Maximum (Minimum) in einer Umgebung des Hopf-Bifurkationspunkts 6 verändert.
In der zweidimensionalen graphischen Darstellung der Realteilverlaufsdarstellung 8 ist auf der horizontalen Achse logk mit einem Wertebereich von "15-30" und auf der vertikalen Achse der Realteil von
g(ik)) mit einem Wertebereich von "0-8" aufgetragen.
In 4 sind der erste Realteilverlauf 9 durch eine durchgezogene Linie, der zweite Realteilverlauf 10 durch eine ge punktete Linie und der dritte Realteilverlauf 11 durch eine fein gepunktete Linie dargestellt.
Der erste Realteilverlauf 9 weist für 15 < logk < 22 einen annähernd konstanten Wert für den Realteil
g(ik)) von 4,6 auf. Für 22 < logk < 23,5 weist der erste Realteilverlauf 9 einen starken Ausschlag nach oben auf. Das Betragsmaximum des ersten Realteilverlaufs 9 ist im Punkt logk=22,7,
g(ik))=7,1 gegeben.
Der zweite Realteilverlauf 10 weist im Bereich 15<logk <22 einen konstanten Wert
g(ik)) = 4,7 auf . Der zweite Realteilverlauf 10 weist im Bereich 22 < logk < 23,5 einen starken Ausschlag nach unten auf und erreicht sein lokales Minimum im Punkt log k = 22,7,
g(ik)) = 1,4.
Der dritte Realteilverlauf 11 verläuft im Bereich 15<logk <22 konstant auf dem Wert
g(ik)) = 4,5. Der dritte Realteilverlauf 11 weist im Bereich 22 < logk < 23,5 ein lokales Minimum bei dem wert log k = 22,4,
g(ik)) = 3,5 auf.
Ab dem Wert logk = 23,5 sinken der erste Realteilverlauf 9, der zweite Realteilverlauf 10 und der dritte Realteilverlauf 11 kontinuierlich bis zu dem Punkt logk = 30,
g(ik)) = 0,2 ab.
Das lokale Maximum des ersten Realteilverlaufs 9, das lokale Minimum des zweiten Realteilverlaufs 10 und des dritten Realteilverlaufs 11 im Bereich 22 < logk < 23,5 weisen auf das Vorhandensein eines Hopf-Bifurkationspunkts 6 hin. Offensichtlich ist der Hopf-Bifurkationspunkt 6 zwischen dem ersten Realteilverlauf 9 und dem zweiten Realteilverlauf 10 bzw. zwischen dem schen dem ersten Fortsetzungsparameter λ₁ und dem zweiten Fortsetzungsparameter λ₂ angeordnet.
Betrachtet man das Monotonieverhalten des Realteils
g(ik)), so stellt man fest, daß bei dem ersten Realteilverlauf 9 vor dem lokalen Maximum eine positive Monotonie, danach eine negative gegeben ist. Bei dem zweiten Realteilverlauf 10 und bei dem dritten Realteilverlauf 11 ist vor dem lokalen Minimum eine negative Monotonie, danach eine positive Monotonie gegeben.
Im Rahmen dieses Dokuments wurden folgende Veröffentlichungen zitiert:

[1]Verfahren und Vorrichtung zur Bestimmung einer periodischen Zustandsbeschreibung eines technischen Systems, welches Schwingungen unterliegt, durch einen Rechner sowie deren Verwendung, 1996, DE 196 40 583 C1 .
[2]R. Neubert, P. Selting, Q. Zheng, Analysis of autonomous oscillators: a multistage approach, in Mathematical theory of networks and systems. Proceedings of the MTNS-98 symposium, held in Padova, Italy, A. Beghi, L. Finesso, and G. Picci, eds., Padova, July 1998, Il poligrafo, pp. 1055-1058.
[3]Q. Zheng, R. Neubert, Computation of periodic solutions of differential algebraic equations in the neighborhood of Hopf bifurcation points, Int. J. of Bifurcation and Chaos, 7 (1997), pp. 2773-2788.
[4]I.N. Bronstein, K.A. Semendjajew, Taschenbuch der Mathematik, 1985, pp. 740-743
[5]I.N. Bronstein, K.A. Semendjajew, Taschenbuch der Mathematik, 1985, pp. 150-159.
[6]M. Tiebout, A Fully Integrated 1.3GHz VCO for GSM in 0.25μm Standard CMOS with a Phasenoise of –142dBc/Hz at 3 MHz Offset, Proceedings 30^th European Microwave Conference, Paris, 2000.
[7]E. Hopf, Abzweigung einer periodischen Lösung von einer stationären Lösung eines Differentialsystems, Bericht der Math.-Phys. Klasse der Sächsischen Akademie der Wissenschaften zu Leipzig, 94 (1942).
[8]G.H. Golub, C.F. van Loan, Matrix Computations, John Hopkins University Press, Baltimore, 1989.
[9]W. Govaerts, A. Spence, Detection of Hopf points by counting sectors in the complex plane, Numer. Math., 75 (1996), pp. 43-58.

1: Flußdiagramm
2: Spannungsgesteuerter Oszillator
3: Bifurkationsdiagramm
4: stabile stationäre Lösungen
5: instabile Lösungen
6: Hopf-Bifurkationspunkt
7: periodische Lösungen
8: Realteilverlaufsdarstellung
9: erster Realteilverlauf
10: zweiter Realteilverlauf
11: dritter Realteilverlauf

Claims

Verfahren zur Bestimmung von Hopf-Bifurkationspunkten (6) einer periodischen Zustandsbeschreibung eines Schwingungen aufweisenden technischen Systems (2), das Gleichungen der folgenden Form genügt: f(x'(t),x(r),λ) = 0 dadurch gekennzeichnet, daß die folgenden Schritte ausgeführt werden: a) Bestimmen eines Werts für den Parameter λ und Berechnen einer stationären Lösung (x,λ) durch Lösen von F(x,λ) = 0,b) Berechnen der Jacobi-Matrizen C, G, c) Aufstellen der Funktion g ~(z) = cT(G + zC)–1b,wobei b, c orthonormierte Zufallsvektoren sind, d) Abspalten des konstanten Terms g_∞ g ~(z) = g∞ + g(z),e) Berechnen der Windungszahl WZ der Funktion g, ausgewertet entlang der imaginären Achse:
f) Bestimmen des Monotonieverhaltens des Realteils
g(ik)) der Funktion g(ik), g) Feststellen, ob ein Hopf-Bifurkationspunkt (6) gefunden wurde, unter Verwendung der Informationen aus den Windungszahlen WZ von g(ik) sowie unter Verwendung der Informationen aus dem Monotonieverhalten von
g(ik)): Falls Hopf-Bifurkationspunkt (6) gefunden wurde: Bestimmen einer Näherung für die Frequenz der Oszillation und gegebenenfalls Durchführen einer inversen Iteration zur genauen Bestimmung der tatsächlichen Frequenz (6), falls kein Hopf-Bifurkationspunkt (6) gefunden wurde: Überprüfen ob λ ≤ λ_max, h) Falls kein Hopf-Bifurkationspunkt (6) gefunden wurde und falls λ ≤ λ_max , Wiederholen der Schritte a) bis f)
Verfahren nach Anspruch 1 dadurch gekennzeichnet, daß als Näherung für die Frequenz einer Oszillation in einer Umgebung des Hopf-Bifurkationspunktes (6) die Stelle k_m gewählt wird, an welcher der Realteil
g(ik)) der Funktion g(ik) seinen steilsten Gradienten hat und wobei die Steigung
der Sekanten von g(ik) durch die Punkte (km, g(ikm)),km–1,g(ikm–1)) auf ihr Betragsmaximum untersucht wird.
Verfahren nach Anspruch 1 oder 2, dadurch gekennzeichnet, daß falls in Schritt g) keine Aussage über das Vorhandensein eines Hopf-Bifurkationspunkts (6) treffbar ist, die Schritte a) bis f) wiederholt werden, wobei in Schritt a) für den neuen Parameter λ das aus dem letzten und aus dem vorletzten Wert des Parameters λ gebildete arithmetische Mittel gewählt wird.
Verfahren nach einem der vorhergehenden Ansprüche, dadurch gekennzeichnet, daß das technische System wenigstens eine elektrische Schaltung aufweist.
Verwendung des Verfahrens nach einem der Ansprüche 1 bis 4 zur elektrischen Schaltungssimulation.
Computerprogrammprodukt sowie Computerprogramm zur Ausführung eines Verfahrens zur Bestimmung von Hopf-Bifurkationspunkten (6) einer periodischen Zustandsbeschreibung eines Schwingungen aufweisenden technischen Systems (2), das so ausgebildet ist, daß ein Verfahren gemäß einem der Ansprüche 1 bis 4 ausführbar ist.
Computerprogramm nach Anspruch 6, das auf einem Speichermedium enthalten ist.
Computerprogramm nach Anspruch 6, das in einem Computerspeicher abgelegt ist.
Computerprogramm nach Anspruch 6, das in einem Direktzugriffsspeicher enthalten ist.
Computerprogramm nach Anspruch 6, das auf einem elektrischen Trägersignal übertragen wird.
Datenträger mit einem Computerprogrammprodukt bzw. Computerprogramm nach Anspruch 6.
Verfahren zum Herunterladen eines Computerprogramms nach Anspruch 6 aus einem elektronischen Datennetz wie beispielsweise aus dem Internet auf einen an das Datennetz angeschlossenen Computer.