Beschreibung
Verfahren und Vorrichtung zur Ermittlung einer Störung eines technischen Systems
Die numerische Simulation elektrischer Schaltungen hat in den letzten Jahren große Bedeutung bei der Entwicklung von Computerchips erlangt. Aufgrund der hohen Kosten für eine Musteranfertigung eines Chips und eines möglichen Redesigns sind Simulatoren unabdingbar geworden. Mit ihnen ist es möglich, prädiktive Aussagen über Betriebsverhalten und Effizienz der modellierten Schaltung an einem Rechner zu erhalten. Nach erfolgreichen Simulationsergebnissen wird üblicherweise ein Chip in Silizium gebrannt.
Unter Verwendung des sogenannten Netzwerkansatzes wird eine Schaltung durch ihre topologischen Eigenschaften, die charakteristischen Gleichungen von Schaltelementen und die Kirch- hoffschen Regeln beschrieben.
Zur Analyse einer Schaltung wird die aus [1] bekannte sogenannte modifizierte Knotenanalyse eingesetzt. Diese führt auf ein differential-algebraisches Gleichungssystem der Form
C(x(t)) • x(t) + f(x(t)) + s(t) = 0. (1)
Allgemein ist unter einem differential-algebraischen Gleichungssystem ein Gleichungssystem des Typs
F x(t), x(t), t
zu verstehen mit einer singulären Jacobi-Matrix F. der par- x tiellen Ableitungen von F nach x(t) .
Das differential-algebraische Gleichungssystem (1) wird auch als quasilineares-implizites Gleichungssystem bezeichnet. Mit x(t) wird ein von einer Zeit t abhängiger Verlauf von Knotenspannungen bezeichnet, mit x(t) dessen Ableitung nach der Zeit t. Mit f(x(t)) wird eine erste vorgegebene Funktion bezeichnet, welche Leitwerte und nicht-lineare Elemente enthält, mit c(x(t)) eine vorgegebene Kapazitätsmatrix und mit s(t) eine zweite vorgegebene Funktion, die unabhängige Spannungsquellen und Stromquellen enthält.
Die Gleichung (1) beschreibt nur den Idealfall für eine Schaltung. In der Praxis läßt sich jedoch Rauschen, d.h. eine Störung der Schaltung, nicht vermeiden. Unter Rauschen versteht man eine ungewollte Signalstörung, die beispielsweise durch thermische Effekte oder die diskrete Struktur der Elementarladung verursacht wird. Aufgrund der fortschreitenden Integrationsdichte integrierter Schaltungen wächst die Bedeutung der prädiktiven Analyse solcher Effekte (Rauschsimulation) .
Bei der Analyse einer Schaltung unter Berücksichtigung von Rauschen kann die Gleichung (1) nodelliert werden als:
C(x(t)) ■ x(t) + f(x(t)) + s(t) + B(t, x(t)) • v(ω, t) = 0. (3)
Hierbei bezeichnet v(ω, t) einen m-dimensionalen Vektor, dessen unabhängige Komponenten verallgemeinertes weißes Rauschen sind, wobei m die Zahl der Störstromquellen angibt. Eine auch als Intensitätsmatrix bezeichnete Matrix ß(t, x(t)) hat die Di- mension n x m. Tritt nur thermisches Rauschen in linearen Widerständen auf, so ist sie konstant.
Wird die Schaltung mit rein linearen Elementen modelliert, so wird die Vorschrift (3) zu einem gestörten linear-impliziten differential-algebraischen Gleichungssystem der Form
C • x(t) + G • x(t) + s(t) + B(t, x(t)) • v(ω, t) = 0. (4)
Unter dem Begriff Index ist im weiteren ein Maß zu verstehen, wie weit sich ein differential-algebraisches Gleichungssystem von einem expliziten gewöhnlichen Differentialgleichungssy- stem "unterscheidet", wie viele Ableitungsschritte erforderlich sind, um aus dem differential-algebraischen Gleichungssystem ein explizites gewöhnliches Differentialgleichungssystem zu erhalten.
Ohne Einschränkung der Allgemeingültigkeit und bei Existenz verschiedener Begriffsdefinitionen des Begriffs "Index", wird im weiteren folgende Definition für den Index eines differen- tial-algebraischen Gleichungssystems verwendet:
Gegeben sei ein differential-algebraisches Gleichungssystem des Typs
x(t), x(t), t = 0 (2:
Existiert eine kleinste natürliche Zahl i, so daß die Gleichungen
x(t), x(t), t 0,
dF x(t), x(t), t
= 0, dt
d^F x(t), x(t), t
= 0 dt (i)
algebraisch in ein System expliziter gewöhnlicher Differentialgleichungen umgeformt werden können, so wird i als der Index des differential-algebraischen Gleichungssystems bezeich- net. Die Funktion F wird dabei als ausreichend oft differenzierbar vorausgesetzt.
Unter einem stochastischen Differentialgleichungssystem ist im weiteren allgemein folgendes Differentialgleichungssystem zu verstehen:
Gegeben sei ein Wiener-Hopf-Prozeß j ^; t e 9ϊc auf einem
Wahrscheinlichkeitsraum (Ω, A, P) zusammen mit seiner kanonischen Filtration {cs; s e [a, b]j . Seien ferner h und G: [a, b] x 9? → 9. zwei NBΓ^I X Bl - Bj -meßbare Zufallsvariablen und X:Ω -. S eine (ca - B)-meßbare Funktion. Ein stochasti- sches Differentialgleichungssystem ist gegeben durch das Itδ- Differential
Xs:Ω→ 3i,
ω i→ X(ω) + J h(u, Xu(ω))dλ(u) + j G(U, Xu)dWu (ω)
[a,s] a
oder symbolisch
dXs = h(s, Xs)ds + G(s, Xs)dWs, s e [a, b]; Xa = X . ( 9 )
Folgendes Verfahren zur Rauschsimulation ist aus [2] bekannt.
Die Vorschrift (4) kann für den Fall einer rein additiven
Störung entkoppelt werden in einen differentiellen und einen algebraischen Teil.
Für den Fall einer rein additiven Störung gilt:
B(t, x(t)) ≡ B(t), do;
d.h. die Intensitätsmatrix hängt nur von der Zeit t ab.
Damit wird (4) zu
C • x(t) + G • x(t) + s(t) + B(t) • v(ω, t) = 0. (iι;
Unter den Annahmen, daß konsistente Anfangswerte
Xdetlto)' Xdetlto) zu einer Anfangszeit t(j gegeben sind und
bei Regularität des Matrixbüschels {λC + G; λ e C}, existiert eine eindeutige Lösung x^e zu (1) in der Form
C • xdet+ G • x^et + s = 0. (12)
Ein Matrixbüschel {λC + G; λ s C} ist regulär, falls ein λg aus C existiert derart, daß gilt:
det(λ0C + G) ≠ 0. (13)
λ
Konsistente Anfangswerte xdet(tθ)' xdet(t0) können dadurch
J gewonnen werden, daß man einen DC-Arbeitspunkt des Systems, welches durch (12) beschrieben wird, bestimmt, also xdet = 0 setzt. Aus [3] ist ein weiteres Verfahren zur Ermittlung kon- f . \ sistenter Anfangswerte xdet(tθ)' xdet(to) bekannt .
Nach einer Transformation, die im folgenden beschrieben wird, gelangt man zu Vorschriften, die zur Vorschrift (11) äquiva- lent sind, der folgenden Form:
y + Fl m + F2 • yM + σM + (B(t) • v(ω, t) - = '14'
und
*3 .H + F4 + σM + (B(t) • v(ω, t))W = .15)
mit transformierten Anfangsbedingungen
'Wdet(to): = -1 .. ,_. * [1]J ( "1 ■ xdet(tθ) :i6:
■[1] [1]
Y det(to): = T • xdet(to) 17;
und
[2] y[2]det(to): = (T X • Xdetfo)) >
•[2] r ] y det(to = T l • xdet(to) :i9.
und mit
B(t): = S • B(t) (20;
Mit Fj_ , i = 1, 2, 3, 4, wird jeweils eine vorgebbare Matrix bezeichnet.
Die Transformation hat zum Ziel, die Vorschrift (11) in ein semi-explizites Differentialgleichungssystem des Typs der
Vorschriften (14) und (15) in der Variable y = f ylAJ, yL2Jj mit
geeigneten Matrizen FJL und einer Funktion σ = σl -, σ- - zu überführen.
Zur Matrix C aus (11) findet man dazu zwei reguläre Matrizen S und T derart, daß die Vorschrift
S • C • T = blockdiag(l, N) (21)
erfüllt ist, wobei mit I eine Einheitsmatrix der Dimension r und mit N eine Nullmatrix der Dimension (n-r) bezeichnet wird.
Mittels Gauß Elimination mit vollständiger Pivotstrategie werden zwei reguläre Matrizen Pj_ und Qi ermittelt mit
PL • C • Qi = CX , (22;
wobei eine Matrix C eine rechte obere Dreicksmatrix ist, deren r erste Diagonalelemente ungleich dem Wert Null sind. Die Matrix C hat ab der (r+l)-ten Zeile einschließlich nur Einträge mit dem Wert Null. Die Matrix Q]_ wird als orthogonale Spalten-Permutationsmatrix gewählt. Die Matrix l ist das Produkt aus einer linken unteren Dreiecksmatrix und einer orthogonalen Zeilen-Permutationsmatrix .
Durch eine Multiplikation mit einer regulären rechten oberen Dreiecksmatrix M von rechts werden alle Einträge der Matrix C]_ oberhalb ihrer Diagonalen eliminiert:
C2: = Ci • M = diag K\ , ... , λm, 0, ... ,0 !23!
(n - r)-mal,
Durch eine Multiplikation mit einer regulären Diagonalmatrix M von links werden alle nichtverschwindenden Diagonalelemente der Matrix C2 auf den Wert 1 transformiert:
C3: = M • C = diag 1, ... ,1, 0, ... ,0 .24 : ^r-mal (n -r)-mal,
Die Matrizen S: = M • P]_ und T: = Qi • M^ leisten das Gewünsch- te.
Durch Setzen von y: = T-1 • x , E: = C3 = S • C • T , F: = S • G • T und σ = S • s in Vorschrift (11) ergibt sich
y+F- y + σ + B- v = 0 (25)
Um die spezielle Struktur der Matrix C3 auszunutzen, wird y in einen ersten Vektor y*- , der die ersten r Komponenten enthält, und in einen zweiten Vektor yfeLi ' , der die restlichen
(n-r) Einträge enthält, unterteilt:
= (y(ll< y 2] [261
Die Matrix F wird in 4 Untermatrizen Fj_ , i = 1, 2, 3, 4 der Dimensionen r x r, r x (n-r) , (n-r) x r, (n-r) x (n-r) aufgeteilt:
JE^
F = : 27 ) \v3fc4j
Für die Matrix E wird eine entsprechende Aufteilung gewählt,
Die Matrix E ist eine Einheitsmatrix der Dimension r x r und die Matrizen E2, E3 und E4 sind Nullmatrizen. Somit zerfällt (25) in die beiden Vorschriften (14) und (15) .
9 Ist die Matrix F4 invertierbar, was genau dann der Fall ist, wenn das System aus Vorschrift (12) den Index 1 besitzt, so kann die Vorschrift (15),
F3 • yW + F4 • yl2] + σt2] + (§(t) • v(ω, t) 2- = 0 ( 15 )
nach yl J aufgelöst werden, was zu folgender Vorschrift führt :
yM = -F4-1 • |F3 • yW + J21 + (§(t) • v(ω, t))[2-j . (28)
Im folgenden werden folgende abkürzende Bezeichnungen eingeführt :
F: = Fλ - F2 • F4 _1 • F3 ; (29 )
s: = σl1-! - F2 • F4 _1 • σl2J ; ( 30 )
B(t, x): = (lr,-F2 • F4 _1) • B(t) . ( 31 )
Mit Ir wird eine Einheitsmatrix der Dimension r, also des
Rangs der Matrix C bezeichnet.
Durch Einsetzen von Vorschrift (28) in Vorschrift (14) ergibt sich:
[1] Ml y + F • yl1J + s + B(t) • v(ω, t) = 0. (32)
Vorschrift (32) kann als ein stochastisches Differentialglei- chungssystem interpretiert werden der folgenden Form:
dY = yW - |F . yW + sldt - B(t)dWt . (33)
10 it Y!. * wird eine Zufallsvariable mit einem Erwartungswert to y t(to) bezeichnet, die eine endliche Varianz aufweist. <Wt; t e 9?Q> ist ein Wiener-Hopf-Prozeß der Dimension der Anzahl der Rauschquellen, allgemein der Anzahl der Störquellen.
Das Verfahren aus [2] geht aus von Vorschrift (33) , deren eindeutiger Lösungsprozeß YfilJ als Itδ-Differential durch die Gleichung
41] - Φt.to 4 )
^ w0 _ J Φt/U • s(u)dλ(u) + J Φt,u • B(u>dw-u : 3
.[to ] t0
mit dem Fundamentalsystem von Lösungen
Φt.tg exp - J Fdλ(u) = exp{(t0 - t)F} :35i [tθt]
gegeben ist.
Bei dem Verfahren aus [2] werden die Erwartungswerte E^ und die zweiten Momente P der Zufallsvariablen Y [ιJl näherungsweise bestimmt. Der Erwartungswert eines Itδ-Integrals ist gleich dem Wert 0. Somit erhält man aus Vorschrift (34) direkt
Et : = E Y w !) = φt,t0 • y^et^o) - I φt,u • έ(u)dλ(u) : 35 j
[tot]
Für alle t löst E^ das gewöhnliche Differentialgleichungssystem
11
xt+ F • x-(- + s = 0; χt0 = &it^o) :36)
Die zweiten Momente
Pt = E 1' ( 37 :
der Zufallsvariablen Y ■■ des Lösungsprozesses genügen für alle t dem Differentialgleichungssystem
= -F • Pt - Pt - s dt Et ST + B(t) • B(t)1 '38'
wobei eine Anfangsbedingung durch
ptn = E' Y .39) tWo • vYtMo
gegeben ist. Bei (38) handelt es sich um ein lineares gewöhnliches Differentialgleichungssystem.
Parallel zur transienten Simulation der Schaltung werden bei dem Verfahren aus [2] durch numerische Integration mit linearen impliziten Mehrschrittverfahren näherungsweise die Erwartungswerte E^ und die zweiten Momente - ermittelt.
Ein Nachteil dieses Verfahrens ist darin zu sehen, daß zur Bestimmung der zweiten Momente P^- für jeden Zeitschritt lineare Gleichungssysteme von quadratischer Ordnung in der Anzahl m der Rauschquellen gelöst werden.
Das Verfahren basiert auf einer manuellen Indexreduktion des differential-algebraischen Gleichungssystems auf ein explizites stochastisches Differentialgleichungssystem, die nicht automatisierbar ist. Außerdem wird nur ein Teil der Rauscheffekte betrachtet. Ferner liefert das Verfahren keine pfadwei-
12 se Information, sondern nur die Momente der Knotenpotentiale, die bei der Indexreduktion erhalten bleiben. Für die Knotenpotentiale, die durch die Indexreduktion unterdrückt werden, liefert dieses Verfahren keine Information.
Ein weiterer Nachteil des Verfahrens aus [2] ist darin zu sehen, daß bei diesem Verfahren die Indexreduktion äußerst ineffizient sowie auch nicht automatisch durchführbar ist, insbesondere da die algebraischen Variablen in dem differen- tial-algebraischen Gleichungssystem nicht berücksichtigt werden. Die Indexreduktion muß bei dem Verfahren aus [2] analytisch manuell durchgeführt werden, da die numerischen Verfahren nicht stabil sind.
Ferner ist es aus [4] bekannt, eine Rauschsimulation einer integrierten Schaltung im Frequenzbereich durchzuführen. Dadurch läßt sich eine Schaltung jedoch nur im Kleinsignalbereich analysieren. Die Voraussetzung eines festen Arbeitspunktes ist jedoch häufig nicht gegeben. Beispielsweise verhindert das Schwingverhalten eines Oszillators in einer Schaltung einen einheitlichen Arbeitpunkt.
In [5] ist eine Schaltungsbeschreibungssprache SPICE beschrieben, mit der eine elektrische Schaltung in einer für einen Rechner verarbeitbarer Form beschreibbar ist.
Aus [6] ist ein Verfahren zur numerischen Behandlung eines stochastischen Differentialgleichungssystems bekannt, die pfadweise Simulation diskreter Approximationen an den Lö- sungsprozeß, das als Runge-Kutta-Schema bezeichnet wird.
Aus [7] ist ein Simulator bekannt, der
• eine Initialisierungseinheit,
• eine Inkrementiereinheit, • eine Einheit zur Aktualisierung eines Schätzwerts, und
• eine Ausgabeeinheit
13 umfaßt. Von der Initialisierungseinheit wird ein Initialisierungswert eines Zustands und einer Funktion, mit denen ein System, welches einer zufälligen Störung unterliegt, beschrieben wird, der Inkrementiereinheit zugeführt. Die Inkre- mentiereinheit verwendet zwei Zufallszahlenfolgen um ein In- krement des Schätzwerts zu bilden, ohne die Funktion selbst zu differenzieren. In der Einheit zur Aktualisierung eines Schätzwerts wird der Schätzwert um das Inkrement erhöht.
Weiterhin sind aus [8] ein Verfahren sowie eine Vorrichtung zur Verbesserung der Genauigkeit eines in einem geschlossenen Regelkreis geregelten Systems bekannt, wobei mindestens zwei stochastische Störsignale beerücksichtigt werden.
Somit liegt der Erfindung das Problem zugrunde, ein Verfahren anzugeben, mit dem die im vorigen beschriebenen Nachteile vermieden werden.
Das Problem wird durch das Verfahren gemäß Patentanspruch 1 und durch die Vorrichtung gemäß Patentanspruch 12 gelöst.
Bei dem Verfahren gemäß Anspruch 1 wird ein technisches System, welches einer Störung unterliegt, mittels eines impliziten stochastischen Differentialgleichungssystems (SDE) be- schrieben. Eine näherungsweise Lösung des Systems wird ermittelt, indem ein diskreter Näherungsprozeß realisiert wird. Der diskrete Näherungsprozeß wird gemäß folgender Vorschrift realisiert:
- (c + hα2G) • XTn + 1 = {- (l - γ)c + h(l + γαi - CX2)G} • XXn +
+{- γC + h(γ - αιγ)GJ • ^^ + +hα2s(τn +ι) + h(l + γoti - α2)s(τn) +
+h(γ - αιγ)s(τn_ι) + +ß(τn) • ΔWn + γB(τn) • ΔWn_ι wobei mit
— C eine erste Matrix,
14
— αi, α , Y vorgegebene Parameter aus dem Intervall [0, 1],
T
— h:= — , eine Schrittweite in einem Ausgangszeitintervall
[0, T] , wobei T ein vorgegebener Wert ist, welches in N Teilintervalle unterteilt ist, — G eine zweite Matrix,
— Xτ - eine Realisierung des Näherungsprozesses an einer
Stützstelle τn+l/
— Xτ eine Realisierung des Näherungsprozesses an einer
Stützstelle τn, — Xτ __. eine Realisierung des Näherungsprozesses an einer
Stützstelle τn-l>
— s(τn+ι) ein erster Wert an der Stützstelle τn +ι,
— s(τn) ein erster Wert an der Stützstelle τn,
— s(τn_ι) ein erster Wert an der Stützstelle τn_ι, — ein Differenzwert ΔWn_ : = Wτ - Wτ »zwischen einem zweiten Wert WT ln an der Stützstelle τ -n l und dem zweiten Wert
Wτ _» an der Stützstelle τn- i
— B(τn) ein zweiter Wert an der Stützstelle τn , bezeichnet wird, und - bei dem durch iterative Lösung des Näherungsprozesses die Störung Xτ - ermittelt wird.
Die Vorrichtung gemäß Patentanspruch 12 weist eine Prozessoreinheit auf, die derart eingerichtet ist, daß ein techni- sches System, welches einer Störung unterliegt, mittels eines impliziten stochastischen Differentialgleichungssystems beschrieben werden. Eine näherungsweise Lösung des Systems wird ermittelt, indem ein diskreter Näherungsprozeß realisiert wird. Der diskrete Näherungsprozeß wird gemäß folgender Vor- schrift realisiert:
15
— (C + hα2G) • XXn + 1 = {- (l - γ)c + h(l + γαι - α2)G} • XXn +
+{- γC + h(γ - αιγ)G} • *τn_1 + +hα2s(τn + ι) + h(l + γαi - α2)s(τn) + +h(γ - αιγ)s(τn_ι) + +ß(τn) • Δwn + γß(τn) • Δwn_x wobei mit
— C eine erste Matrix,
— αi, α2, γ vorgegebene Parameter aus dem Intervall [0, 1],
T — h:= — , eine Schrittweite in einem Ausgangszeitintervall N -
[0, T] , wobei T ein vorgegebener Wert ist, welches in N
Teilintervalle unterteilt ist,
— G eine zweite Matrix,
— Xτ - eine Realisierung des Näherungsprozesses an einer Stützstelle τn+ ,
— Xτ eine Realisierung des Näherungsprozesses an einer
Stützstelle τn ,
— Xτ __. eine Realisierung des Näherungsprozesses an einer
Stützstelle τn_l' — s(τn+l) ein erster Wert an der Stützstelle ιn + ι ,
— s(τn) ein erster Wert an der Stützstelle τn ,
— s(τn-l) ein erster Wert an der Stützstelle τn_ι,
— ein Differenzwert ΔWn_ι: = Wτ - Wτ - zwischen einem zweiten Wert Wτ an der Stützstelle τn und dem zweiten Wert Wtn-1 an der Stützstelle τn_ι
— B(τn) ein zweiter Wert an der Stützstelle τn , bezeichnet wird. Durch iterative Lösung des Näherungsprozesses wird die Störung ^τn+ι ermittelt.
Die Erfindung verwendet direkt die implizite Struktur des technischen Systems, repräsentiert durch ein implizites Differentialgleichungssystem. Durch die Erfindung wird die Ermittlung der Störung erheblich beschleunigt, da die Dünnbe- setztheit der Matrizen C und G ausgenutzt werden kann. Die
16 numerisch instabile und aufwendige Transformation des Differentialgleichungssystems in die entkoppelte Form entfällt.
Durch die Erfindung ist es erstmals möglich, eine Störung auch bei einer singulären Matrix C zu ermitteln.
Vorteilhafte Weiterbildungen der Erfindung ergeben sich aus den abhängigen Ansprüchen.
Die Erfindung kann in verschiedensten Anwendungsgebieten eingesetzt werden, in denen ein technisches System gestört ist und durch ein System differential-algebraischer Gleichungen beschrieben werden kann.
Beispielsweise können Störungen (Rauschen) in einer elektrischen Schaltung ermittelt werden. Die Erfindung eignet sich ferner zur Anwendung in einem mechanischen Mehrkörpersystem oder in einem allgemeinen physikalischen System, einem chemischen System oder auch einem physikalisch-chemischen System, deren jeweilige Modellierung auf ein System differentialalgebraischer Gleichungen führt.
Ohne Einschränkung der Allgemeingültigkeit wird im weiteren die Erfindung anhand eines Ausführungsbeispiels einer Rausch- Simulation einer elektrischen Schaltung beschrieben, das in den Figuren dargestellt ist.
Es zeigen
Figur 1 ein Ablaufdiagramm, in dem die einzelnen Verfahrensschritte dargestellt sind;
Figur 2 eine Vorrichtung, mit der das Verfahren durchgeführt wird; Figur 3 eine Skizze einer elektrischen Schaltung eines Differentiators; Figur 4 eine Darstellung eines Simulationsergebnisses.
17
In Figur 2 ist in Form einer Blockskizze ein Rechner R dargestellt, mit dem das im weiteren beschriebene Verfahren durchgeführt wird. Der Rechner R weist eine Prozessoreinheit P auf, die derart eingerichtet ist, daß die im weiteren be- schriebenen Verfahrensschritte durchgeführt werden können. Die Prozessoreinheit P ist über einen Bus B mit einem Speicher SP verbunden. Der Rechner R ist über eine erste Schnittstelle I/O 1 mit einer Tastatur TA und über eine zweite Schnittstelle I/O 2 mit einer Maus MA verbunden, mit denen jeweils ein Benutzer des Rechners R Eingaben vornehmen kann. Ferner ist der Rechner R mit einem Bildschirm BS verbunden, auf dem dem Benutzer Ergebnisse des Verfahrens dargestellt werden.
In dem Speicher SP ist die zu analysierende Schaltung S
(vgl. Figur 3) in Form einer Schaltungsbeschreibungssprache gespeichert. Als Schaltungsbeschreibungssprache wird das Netzlistenformat von SPICE eingesetzt.
Ein Ersatzschaltbild der Schaltung S für dieses Ausführungsbeispiel ist in Figur 3 dargestellt. Die Schaltung S stellt in ihrer Funktionalität einen Differentiator dar. Eine unabhängige Spannungsquelle V(t) ist zwischen einem ersten Knoten N]_ und einem Massepotential NQ angeordnet. Zwischen dem er- sten Knoten und einem zweiten Knoten N2 ist eine Kapazität CQ vorgesehen. Der zweite Knoten N ist über einen Widerstand R mit einem dritten Knoten N3 verbunden. Ferner ist zur Modellierung eines thermischen Rauschens in dem Widerstand R eine zu dem Widerstand R parallel geschaltete Rausch- Stromquelle RS mit der Stromstärke ΔlR vorgesehen. Ein erster Eingang El eines Operationsverstärkers OV ist mit dem zweiten Knoten N2 verbunden, ein zweiter Eingang E2 des Operationsverstärkers OV mit dem Massepotential NQ . Ein Ausgang A des Operationsverstärkers OV ist mit dem dritten Knoten N3 verbunden.
18
Die Schaltung S wird in einem ersten Schritt 101 in Form des Netzlistenformats von SPICE in dem Speicher SP gespeichert (vgl. Figur 1) .
In einem zweiten Schritt 102 wird für die Schaltung S eine modifizierte Knotenanalyse durchgeführt. Das Ergebnis der modifizierten Knotenanalyse ist das der Schaltung S zugehörige Gleichungssystem.
In Matrixschreibweise ergibt sich für die Schaltung S das linear-implizite Differentialgleichungssystem des Indexes 1 mit rein additiver Störung
( ■ λ ul f c0 - c0 0 0 0
- c0 c0 0 0 0 u2
0 0 0 0 0 +
U3
0 0 0 0 0
. o 0 0 0 o , ~~Q
0 0 0 1 0
1 1 ( Ul ^ 4kT
0 0 0 Δf f 0 R R u2 0
1 1 4kT
+ 0 0 - 1 U3 + + Δf v(t, ω) 0 R R R 0
0 A 1 0 0
V-'-ampl/ 0 v(t)J
1 0 0 0 0
Es werden mit
- u, u2, U3 die Spannungen an den Knoten N_, N2, N3 ,
- IQ der durch den ersten Knoten N]_ und durch die
Spannungsquelle V(t) fließende Strom,
- Ia pl der durch den dritten Knoten N3 und den
Operationsverstärker OV fließende Strom,
- R der Wert des Widerstand R,
- CQ der Wert der Kapazität CQ ,
19
- A ein Verstärkungsfaktor des Operationsverstärkers OV, bezeichnet .
In dem Modell des thermischen Rauschens ist die Stromstärke ΔlR des Stroms, der durch die Rauschstromquelle RS fließt, durch die Gleichung
4kT
ΔIR Δf • v(ω, t) 40;
R
angenommen.
Hierbei beschreibt v(ω, t) als weißes Rauschen den treibenden verallgemeinerten stochastischen Prozeß, also die Störung des technischen Systems. Die Größen k, T und Δf sind als Kon- stanten vorgegeben.
Allgemein ergibt sich also ein nichtentkoppeltes linearimplizites Differentialgleichungssystem mit rein additiver stochastischen Störung (Rauschen) der folgenden Form:
C • x(t) + G • x(t) + s(t) + ß(t) • v(ω, t) = 0. (4i:
Zur Herleitung des Verfahrens, nicht aber für das Verfahren selbst, wird angenommen , daß die Kapazitätsmatrix C inver- tierbar ist.
Durch Multiplikation der Gleichung (41) mit der inversen Ka- pazitätsmatrix C erhält man
x(t) = -C_1 • {G • x(t) + s(t)} - C_1 • B(t, x(t)) • v(ω, t) . (42)
Die Vorschrift (42) wird als stochastisches Differentialgleichungssystem der Form
dXt = Xt0 - C_1 • {G • Xt + s}dt - C_1 • B(t)dWt (43)
20
interpretiert, wobei X^ eine Zufallsvariable mit Erwartungswert det(tθ) bezeichnet.
Auf das Differentialgleichungssystem (37) wird folgendes Verfahren, das in [6] beschrieben ist, angewendet (Das Verfahren wird als implizites starkes Zweischrittverfahren der Ordnung 1 für ein stochastisches Differentialgleichungssystem bezeichnet) :
Ein numerisches Verfahren liefert für jede natürliche Zahl N Realisierungen eines diskreten Näherungsprozesses |xs; s = XQ, ...,τn|, indem es für gegebene ω aus Ω an den jeweiligen Stützstellen τ Approximationen x(ω, τj_) für die Wer- te x(ω, τ berechnet.
Um einen stetigen Näherungsprozeß xs; s e [0, T]| zu erhalten, interpoliert man die erhaltenen Werte z.B. stückweise linear. Die sich ergebenden Pfade können anschließend mit statisti- sehen Methoden analysiert werden.
In [6] ist folgende Familie von impliziten starken Zweischrittverfahren der Ordnung 1 zur Bestimmung der Komponente k der Realisierungen |xτ . ; i = 0, ...,N> angegeben:
*n+l = (l " Yk)Xn + YkXΪ^
α2,kfk(τn+l' Hi+l) + i1 + γkαl,k α2,k + h , . y : 44 :
+ Yk(l - αi,k)fk(τn-l^τn_ι)
+ v + YkV^!
Hierbei steht eine Funktion ohne Argumente für eine Auswertung am Punkt lτn, Xτ ). Die reellen Parameter α _k, α2, und γk sind aus dem Intervall [0, 1] gewählt. Ferner ist die Grö- ße Vn gegeben durch
21
m m = ∑gk,jΔW^ + ∑ l lgk'j2l x (45)
. - . . \Jl'LI2}'τΩ.'τΩ. + l
3 = 1 31,32
wobei
Δ n:= WTn + 1 - WTn (46)
die Zuwächse des Wiener-Prozesses beschreibt und L eine abkürzende Schreibweise für den Operator
d . d I?:= ∑^'^— j = l, ...,m (47) k=l &
ist.
Die Größe b beesscchhr:eibt das mehrfache Itδ- bl>32\ τn'' τn + l
Integral
τn + 1
s2 dW^dW
112 ) hlr32\
ττifi
:n + l
' I f J ^ si ≤2 :48
^n τ
n
h ist gegeben durch die Vorschrift
T h:= -. (49)
N
Wählt man für alle k = 1, ... , n die Parameter α ,k, α2,k uni γk, so vereinfacht sich Vorschrift (44) für den Fall rein additiver Störung in Vektorschreibweise zu
*n + 1 = (X " Y)*n + YXτn_!
+h{α2f(τn + ι, Xτn + 1) + i1 + Yαl " α2)f + γ(l " αl)φn-l> Xτn_ι)} / +g ■ ΔWn + γg • ΔWn_ι
22
:5o;
da die Funktion g nur von der Zeit abhängt. Vorschrift (50: wird auf das Differentialgleichungssystem (43) angewendet.
Es ergibt sich somit:
X = (l - γ)χ + γx τn+l '' τn ' τn-l α2C_1(G • Xχn + 1 + s(τn + 1)) +
-h +(l + γαi - α2)c_1(G • Xtn + s(τn)) .5i: +(γ - αιγ)c_1(G • Xτn_1 + s(τn_j_))
-C-1B(τn)ΔWn - γC_1B(τn)ΔVik-ι
Im nächsten Schritt wird (40) mit -C multipliziert.
Es ergibt sich somit:
- C {*τn + l " t1 " )xτn - YXχ n_ι} = «2{G • Xχn + 1 + s(τn + ι)} + (l + γαi - CX2){G • XXn + s(τn)} P2] + (γ - αιγ){G • XXn_1 + s(τn_ι)} +ß(τn) • Δwn + γß(τn) • Δwn_χ
Durch Zusammenfassen der Terme in Xτ - auf der rechten Seite von (41) ergibt sich:
23 - (c + hα2G) • XXn + 1 - {- (l - γ)c + h(l + γαi - CX2)G} • XTn +
+{- γC + h(γ - αιγ)G} • ^ττχ _ 1 + +hα2s(τn + ι) + h(l + γαi - α2)s(τn) + +h(γ - αιγ)s(τn _ !) + +ß(τn) • ΔWn + γß(τn) • ΔWn _ !
( 53 )
Für den Fall des Indexes 0 ist durch Vorschrift (53) eine Vorschrift angegeben, mit der die Störung des Systems ermittelt werden kann. Vorschrift (53) kann auch im Fall des Indexes 1 angewendet werden.
Durch iterative Lösung des Näherungsprozesses wird in einem weiteren Schritt 103 die Störung Xχn+1 ermittelt. Die Bestimmung der Störung τn+ι wird iterativ in einer Schleife für n = 0,1, ... , N - 1 durchgeführt, wodurch der Verlauf der Störung Xτ - zu den jeweiligen Zeitpunkten τπ ermittelt wird.
Abhängig von der ermittelten Störung wird eine Schaltung S modifiziert, so daß die vorgebbaren Bedingungen z.B. hinsichtlich der Störanfälligkeit der Schaltung S erfüllt sind.
Die modifizierte Schaltung wird in einem letzten Schritt in Silizium gebrannt.
Figur 4 zeigt als Ergebnis des von dem Rechner R durchgeführten Verfahrens, das im weiteren als Programm angegeben ist. Es ist dargestellt ein numerisch simulierter Lösungspfad RLP des Spannungsverlaufs im dritten Knoten N3 unter Berücksichtigung des Rauschens. Ferner ist zum Vergleich der ideale Lösungspfad ILP, d.h. der Lösungspfad ohne Berücksichtigung des Rauschens angegeben. Außerdem ist der der Verlauf V der Ein- gangsgröße V(t) dargestellt.
24 Für das Verfahren wurden folgende Parameter verwendet:
2,5 • 10 -11 [s]
0 für 0 < t < 5 • 10 -9 [V]
2 • 10- t - 1 10" für 5 • 10 -9 < t < 10 • 10" v(t) = lO-2 für 10 • 10 -9 < t < 15 • 10-9
-2 • 106w • t + 4 • 10" für 15 • 10 -9 < t; << 20 • 10"
0 für 20 • 10" < t < 25 • 10,-9
C0 = 1 - 10 -12 [F] R = 1 • 104 [Ω]
Δf = 1 [Hz :
Im weiteren sind einige Varianten bzw. Verallgemeinerungen des oben beschriebenen Ausführungsbeispiels dargestellt:
Die Kapazitätsmatrix kann ohne weiteres auch singulär sein. Somit ist es erstmals möglich, auch eine singuläre Kapazitätsmatrix zu berücksichtigen.
Ferner können auch weitere in [6] beschriebene Verfahren zur pfadweisen Approximation eines stochastischen Differentialgleichungssystems eingesetzt werden, um z.B. eine höhere Konvergenzordnung zu erreichen.
Die Erfindung ist keineswegs auf die Ermittlung thermischen Rauschens beschränkt. Im Gegenteil kann jedes Modell einer Störung im Rahmen der Erfindung berücksichtigt werden, welches durch weißes Rauschen beschrieben oder approximiert werden kann, z.B. das Schrotrauschen in Stromquellen. Bei Berücksichtigung des Schrotrauschens wird die Stromquelle modelliert durch eine parallel geschaltete, zufallsabhängige Stromquelle, deren Stromstärke ΔIQ der Vorschrift
ΔIQ = ^2eIcietΔf ' v(ω' t)
genügt <
25
Im weiteren ist eine Realisierung des Verfahrens zur Rauschsimulation in der Programmiersprache FORTRAN77 angegeben.
PROGRAM integr C
C Verfahren zur numerischen Loesung linear-impliziter
C Differentialgleichungssysteme mit rein additiver Störung.
C
IMPLICIT NONE INTRINSIC dble, real
EXTERNAL setall, gennor, dgesv, giveB, gives EXTERNAL a_smsm, m_mv, v_to_v, m_sv, a_vvvv REAL gennor C C Konstanten- und Variablendeklarationen und Initialisierung C
INTEGER n, , anz, i, j, info, initl, init2 C
C n ist die Dimension des Systems, m die Anzahl der C Rauschquellen
C anz ist die Anzahl der betrachteten Teilintervalle, C i und j sind Hilfsvariablen für Schleifen, C info fragt den Return-Code C initl und init2 sind die Samen fuer den C Pseudozufallszahlengenerator C
PARAMETER (n = 5, m = 1, anz = 1000) PARAMETER (initl = 23556, init2 = 4285979) C DOUBLE PRECISION lrand, rrand, h, tauakt, alpha C
C [lrand, rrand] bezeichnet das betrachtete Intervall C h die Schrittweite, tauakt den aktuellen Zeitpunkt, C alpha den Parameterwert alpha2 aus (53), wobei alphal = 0.0 C und gamma = 0.0 gesetzt wurde. C
PARAMETER (alpha = 0.9)
26
PARAMETER (lrand = O.OdO, rrand = 2.5d-8) C
DOUBLE PRECISION xakt (n) , d(n), salt (n) , sakt (n)
DOUBLE PRECISION C(n,n), G(n,n), B(n,m), A(n,n) DOUBLE PRECISION hl_n (n) , h2_n(n), h3_n(n), h4_n(n)
DOUBLE PRECISION hl_m(m), hl_nn(n, n) C
C xakt beschreibt den Wert des Lösungsprozesses am Zeitpunkt C tauakt, C C, B und G sind die Matrizen des Problems der Vorschrift C (3)
C A und d dienen zum Gleichungsaufbau C salt und sakt beschreiben s(taukt) bzw. s(takt - 1) C h*_* bezeichnen Hilfsvektoren und -matrizen C
C Konsistenten Anfangsvektor setzen xakt(l) = O.OdO xakt(2) = O.OdO xakt (3) = O.OdO xakt (4) = O.OdO xakt (5) = O.OdO C
C Definition der Matrix C (mit Kapazität = 1* 10Λ{-12} C Farad) C(l,l) = 1.0d-12
C(2,l) = -1.0d-12
C(3,l) = O.OdO
C(4,l) = O.OdO
C(5,l) = O.OdO C(l, 2) = -1.0d-12
C(2,2) = 1.0d-12
C(3,2) = O.OdO
C(4,2) = O.OdO
C(5,2) = O.OdO
C(l, 3) = O.OdO
C(2,3) = O.OdO
C(3,3) = O.OdO
27
C(4 .3) =_ O.OdO
C(5,3) _= O.OdO
Cd, 4) = O.OdO
C(2,4) _= O.OdO
C(3,4) = O.OdO
C(4,4) = O.OdO
C(5,4) = O.OdO
Cd, 5) = O.OdO
C(2,5) = O.OdO
C(3,5) = O.OdO
C(4,5) = O.OdO
C(5,5) = O.OdO
C
C Definition der Matrix G (mit Rl = 10000, A = 300})
G(l,l) = O.OdO
G(2d) = O.OdO
G(3 ,1) = O.OdO
G(4 ,1) = O.OdO
G(5 ,1) = l.OdO
G(l ,2) = O.OdO
G(2 , 2 ) = l.Od-4
G(3 , 2 ) = -1.0d-4
G(4 .2) _= 3.0d2
G(5 -2) = O.OdO
G(l -3) = O.OdO
G(2 3) = -1.0d-4
G(3 3) = 1.0d-4
G(4 3) = l.OdO
G(5 3) = O.OdO
G(l, 4) = l.OdO
G(2, 4) = O.OdO
G(3, 4) = O.OdO
G(4, 4) = O.OdO
G(5, 4) = O.OdO
6(1, 5) = O.OdO
G(2, 5) = O.OdO
G(3, 5) = -l.OdO
28
G ( 4 , 5 ) = O . OdO G ( 5 , 5 ) = O . OdO C
C Pseudozufallszahlengenerator initialisieren CALL setall (initl, init2) C
C Schrittweite berechnen C h = (rrand - lrand) / anz C
C Zeitpunkt und x-Wert initialisieren C tauakt = lrand C C Ausgabedatei öffnen
OPEN(9, FILE = OpmitlOOO', FORM = 'FORMATTED') C Anfangszeitpunkt und -wert in die Ausgabedatei schreiben
WRITE(9, 42) tauakt, (xakt(i), i = 1, n) C "alten" s-Wert merken CALL gives(salt, n, tauakt) C
C Äussere Schleife: Entspricht der Behandlung eines C Teilintervalls
DO 10, i = 1, anz C Aufbauen des Systems A*X=b C
C Aufbauen von A C A = -C -alpha * h * G
CALL a_smsm(C, G, -l.OdO, -alpha * h, A, n, n) C B aufbauen
CALL giveB(B, n, m, tauakt) C Zeitpunkt erhöhen tauakt = tauakt + h C sakt aufbauen — salt liegt vor Schleifeneintritt vor CALL gives(sakt, n, tauakt)
C "alten" s-Wert merken
CALL v to v(sakt, salt, n)
29
C d aufbauen
C
C hl_nn = -C + (1- alpha) * h * G
CALL a_smsm(C, G, -l.OdO, (1-alpha) * h, hl_nn, n, n) C hl_n = hll_nn * xakt
CALL m_mv(hl_nn, xakt, hl_n, n, n) C h2_n = alpha * h * sakt
CALL m_sv(sakt, alpha*h, h2_n, n) C h3_n = (1 - alpha) * h * salt CALL m_sv(salt, (1- alpha) * h, h3_n, n)
C hl_m = DeltaW
DO 20, j = 1, m hl_m(j) = DBLE (gennor (0.0, SQRT (REAL (h) ) ) ) 20 CONTINUE C h4_n = B * hl_m
CALL m_mv(B, hl_m, h4_n, n, m) C d = hl_n + h2_n + h3_n + h4_n
CALL a_vvvv(hl_n, h2_n, h3_n, h4_n, d, n) C Aufrufen des Gleichungsl'Osers (hl_n ist nur Dummy) CALL dgesv(n, 1, A, n, hl_n, d, n, info)
C Neuen Wert von xakt setzen
CALL v_to_v(d, xakt, n) C aktuellen Zeitpunkt und x-Wert in die Ausgabedatei schrei C ben WRITE(9, 42) tauakt, (xakt(j), j = 1, 3)
42 FORMAT(E16.6,E16.6,E16.6,E16.6,E16.6, E16.6) 10 CONTINUE C
C Ausgabedatei schliessen CLOSE(9, STATUS = 'keep') STOP END C
C
SUBROUTINE giveB(outB, n, m, tauakt) C
30
C Liefert in outB den Wert der (n x m) -Matrix B
C aus Vorschrift (3) zum Zeitpunkt tauakt
C
IMPLICIT NONE
INTRINSIC dsqrt C C Dummy Argumente
INTEGER n,m
DOUBLE PRECISION outB(n, m) , tauakt
C Lokale Variable
DOUBLE PRECISION k, T, deltaf C k bezeichnet die Boltzmann-Konstante, T die absolute C Temperatur, deltaf die Rauschbandbreite PARAMETER(k = 1.308d-23, T = 300, deltaf = l.OdO)
C Zeitunabhängige Definition der Matrix B aus Beispiel 3.4 C (mit Rl = 3000, R2 = 4000, R3 = 5000 Ohm C T = 300 K, k = 1.3807 * 10Λ{-23} jK Λ{-1}, f = 10) C outB(l,l) = O.OdO outB(2,l) = dsqrt (4*k*T*deltaf/10000. OdO) outB.3,1) = -dsqrt (4*k*T*deltaf/10000. OdO) outß (4,1) = O.OdO outß(5, 1) = O.OdO RETURN
END
C SUBROUTINE gives (outs, n, tauakt) C
C Liefert in outs den Wert des (n) -Vektors s C aus Vorschrift (3) zum Zeitpunkt tauakt C C
IMPLICIT NONE INTRINSIC sin
31
C
C Dummy Argumente INTEGER n
DOUBLE PRECISION outs (n) , tauakt C
C Lokale Variable
INTEGER i C i beschreibt die aktuelle Zeilenposition von outs DO 10, i = l, n - 1 outs(i) = O.OdO
10 CONTINUE
IF (O.OdO .le. tauakt .and. tauakt .le. 5.0d-9) THEN outs(n) = O.OdO ELSE IF (5.0d-9 .lt. tauakt .and. tauakt .le. lO.Od-9; THEN outs(n) = -(2.0d6 * tauakt - 1.0d-2) ELSE IF (10.0d-9 .lt. tauakt .and. tauakt .le. 15.0d-9) THEN outs(n) = -(10.0d-3) ELSE IF (15.0d-9 .lt. tauakt .and. tauakt .le. 20.0d-9; THEN outs(n) = -(-2.0d6 * tauakt + 4.0d-2) ELSE IF (20.0d-9 .lt. tauakt .and. tauakt .le. 25.0d-9) THEN outs(n) = O.OdO
END IF RETURN END
32
Im Rahmen dieses Dokuments wurden folgende Veröffentlichungen zitiert:
[1] A. F. Schwarz, Computer-Aided design of microelectronic circuits and Systems, vol. 1, Academic Press, London, ISBN 0-12-632431-X, S. 185 - 188, 1987
[2] A. Demir et al, Time-domain non-Monte Carlo noise
Simulation for nonlinear dynamic circuits with arbitrary excitations, IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, Vol. 15, No . 5, S. 493 - 505, May 1996
[3] B.J. Leimkuhler et al, Approximation methods for the consistent initialization of differential-algebraic Systems of equations, SIAM J. Numer. Anal., Vol. 28, S. 205 - 226, 1991
[4] L. 0. Chua und P. M. Lin, Computer aided design of electronic circuits, Prentice Hall, Englewood Cliffs, 1975, ISBN 0-13-165415-2, S. 596
[5] W. Nagel, SPICE2 - a Computer program to simulate semiconductor circuits, Tech. Report, UC Berkeley, Memo ERL-M 520, 1975
[6] P.E. Kloeden und E. Platen, Numerical solution of stochastic differential equations, Springer Verlag, Berlin, New York, ISBN 3-540-54062-8, S. 412, 1992
[7] US 5 646 869
[8] US 5 132 897