Implementierung

Iterative Berechnung der Größen

Wenn die Schätzung online erfolgen soll, dann kommt zu jedem Abtastschritt ein neuer Datensatz hinzu.

Es müssen aber nicht in jedem Abtastschritt Produkte immer größer werdender Matrizen berechnet werden, sondern man kann die Berechnung iterativ ausführen.

Die MatrixΨ und der Vektory werden hier nur in den Kombinationen L_k=Ψ^T_kΨk

und

r_k=Ψ^T_kr_k

verwendet. Diese können iterativ über L_k=

(0 k=0 L_k−1+ψ_kψ^T_k k>0 und

r_k=

(0 k=0 r_k−1+˜y_kψ_k k>0

berechnet werden (siehe auch Anhang B.3).

Zum Bestimmen der Parameter über θˆk=L⁻_k¹r_k

muss also in jedem Schritt (wenn alle Eingänge ausreichend angeregt wurden) eine (p+1)×(p+1) -Matrix invertiert bzw. ein lineares Gleichungssystem(p+1)-ter Ordnung gelöst werden.

Für den Fallp=2ergibt sich der in Abbildung 6.3 dargestellte Ablaufplan. Solange noch keine Änderung der sekundären Eingangsgröße u₂ vorliegt, kann θ₂ auch noch nicht geschätzt werden. Die Schätzung vonθ₁ erfolgt dann über eine Zwei-Parameter-Schätzung. Siehe dazu auch Abschnitt 6.4.4.

6.4 Identifikation während eines Auftrags 101

Neue Messdaten ψ, ˜y(k)

σ²_ˆ

θ1 = σ²_ˆ

θ2 = ∞

L = L+ψψ^T r = r+ψ˜y

rank(L) =3? L_T = L[1 2],[1 2]

r_T=r_{[1 2]}−L_{[1 2],3}·θˆ_2,r

u₂=const∧ rank(L_T) =2?

[∗,θˆ₁,θˆ₂] = L⁻¹r [∗,σ²_ˆ

θ₁,σ²_ˆ

θ₂] = diag(σ²_yL⁻¹) [∗,θˆ₁] = L⁻¹_T r_T [∗,σ²_ˆ

θ₁] = diag(σ²_yL⁻¹_T )

σ²_ˆ

θ1≤σ²_ˆ

θ1,r,max? θˆ_1,r = ˆθ₁

σ²_ˆ

θ2≤σ²_ˆ

θ2,r,max? θˆ_2,r = ˆθ₂

θˆ_0,r=L⁻¹_1,1·(r₁−L_{1,[2 3]}·[ ˆθ_1,r θˆ_2,r]^T)

Reglerparameterθˆ_0,r,θˆ_1,r,θˆ_2,r Ja

Nein

Ja

Nein

Ja

Nein

Ja

Nein

Initialisierung bei neuem Auftrag

L=0 r=0

θˆ_0,r,θˆ_1,r,θˆ_2,r laden

Abbildung 6.3:Ablaufplan Drei-Parameter-Schätzung

Identifikation bezüglich des Nennpunkts

In Abbildung 6.4 ist beispielhaft dargestellt, weshalb es günstiger ist, die Identifikation bezüglich der Nenngrößen, d. h. dem Zusammenhang

y_k0 =θ₀+θ₁·(u_1,k0−u_1,N) +θ₂·(u_2,k0−u_2,N) +· · · durchzuführen, und nicht den Zusammenhang

y_k0 =θ₀₀+θ₁·u_1,k0+θ₂·u_2,k0+· · · zugrunde zu legen.

Dazu ist in Abbildung 6.4 für den Fallp=1die Schätzung der Geraden mit zwei Messungen beiu₁=40 und 60 dargestellt. Die durchgezogene Linie entspricht dem tatsächlichen Zusammenhang. Durch das Messrauschen entstehen bei den beiden Eingangswerten von 40 und 60 Unsicherheiten, die durch die kleinen senkrechten Striche dargestellt sind. Der Bereich, in dem die damit geschätzte Gerade liegt, ist durch den grauen Bereich angedeutet. Es ist zu erkennen, dass dieser Bereich beiu₁=50deutlich kleiner ist als beiu₁=0. D. h. die Varianz des Schätzwertes{θˆ₀₀}rwäre deutlich größer als die des Schätzwertes {θˆ₀}r. Da man sich in der Anwendung immer um den Nennpunkt aufhalten wird, ist die Varianz von {θˆ₀}r aussagekräftiger. (Vgl. dazu auch die Bemerkung von [F^REUNDet al., 2006, S. 42, 50] bezüglich der Interpretierbarkeit vonθˆ₀.)

{θˆ₀₀}r

{θˆ₀}r

0 10 20 30 40 50 60 70

10 20 30 40 50 60

u₁ y

Abbildung 6.4:Vergleich vonθ₀undθ₀₀

Skalierung der Eingangsgrößen

Um Einträge gleicher Größenordnung in derL-Matrix und damit günstigere Bedingungen für die nume-rische Berechnung zu erhalten, werden die Eingangsgrößen normiert. Dazu wird die Gleichung

y_k0 =θ₀+θ₁·(u_1,k0−u_1,N) +θ₂·(u_2,k0−u_2,N) in

y_k0 =s₀·θ₀

s₀ + (u_1,k0−u_1,N)·s₁·θ₁

s₁ + (u_2,k0−u_2,N)·s₂·θ₂ s₂ umgeschrieben. Nun werden die Größen ^θ_s⁰

0, ^θ_s¹

1 und ^θ_s²

2 als neue Parameter θ₀⁰, θ₁⁰ und θ₂⁰ aufgefasst.

Damit wird der Parametervektor θ⁰=”

θ₀⁰ θ₁⁰ θ₂⁰—_T

gesucht, wobei der Eingangsvektor ψ⁰_k^T0 =”

s₀ (u_1,k0−u_1,N)·s₁ (u_2,k0−u_2,N)·s₂— lautet.

6.4 Identifikation während eines Auftrags 103

Algorithmen zur numerischen Lösung

Im Ablaufplan in Abbildung 6.3 ist sowie der Gleichung (6.4) bzw. (6.3) wird der Schätzwert über die Inverse der MatrixL=Ψ^TΨ bestimmt,θˆ =L⁻¹r. Dies ist jedoch nicht „wörtlich“ zu verstehen, sondern der Vektor der Schätzwerte wird numerisch über die Lösung des Gleichungssystems

Lθˆ =r,

der sogenannten Normalgleichung, bestimmt. Diese Lösung wird, da L, wenn Ψ vollen Spalten-rang hat, symmetrisch und positiv definit ist, über eine Cholesky-Zerlegung von L berechnet.

[DAHMENund REUSKEN, 2008, S. 127] Dabei kann die Inverse von L, die zur Bestimmung der Varian-zen der Schätzwerte benötigt wird, parallel zur Lösung θˆ berechnet werden, d. h. es muss nur einmal die Cholesky-Zerlegung durchgeführt werden. [D^AHMENund R^EUSKEN, 2008, S. 80]

Grundsätzlich ist bekannt, dass die Bestimmung der Lösung des LS-Problems über die Berechnung der Lösung vonLθˆ =r bzw.Ψ^TΨθˆ =Ψ^T˜y numerisch ungünstig sein kann und bei schlecht konditionierten Matrizen Lzu Problemen führt, siehe z. B. [D^AHMENund R^EUSKEN, 2008, S. 128]. Darüber hinaus wird an der genannten Stelle auch darauf hingewiesen, dass schon die Bestimmung der MatrixL=Ψ^TΨ für großekaufgrund von Rundungsfehlern kritisch sein kann.

Prinzipiell wird in der Literatur häufig empfohlen, die Lösung des LS-Problem direkt, basierend auf einer QR-Zerlegung der Eingangsmatrix Ψ zu bestimmen, siehe z. B. [DAHMENund REUSKEN, 2008, S. 129], [VONFINCKENSTEINet al., 2002, S. 222].

Für den vorliegenden Zweck ist die Bestimmung der Lösung des LS-Problems über die Normalgleichung jedoch unkritisch. Dies liegt daran, dass ein schlecht konditioniertesLdarauf hinweist, dass zumindest ein Parameter nur sehr schlecht, d. h. mit hoher Varianz geschätzt werden kann. Es besteht damit keine Notwendigkeit, die Schätzwerte bei einer schlecht konditionierten MatrixLzu bestimmen, da diese dann ohnehin nicht berücksichtigt werden dürfen. Es muss also nur sichergestellt sein, dass die schlechte Kondition vonLerkannt wird, bevor die Schätzwerte bestimmt werden.

Der Rang einer Matrix kann numerisch über die Berechnung der Singulärwerte bestimmt werden, siehe z. B. [D^AHMENund R^EUSKEN, 2008, S. 150]. Die berechneten Singulärwerte können hierbei auch genutzt werden, bei zu schlecht konditionierten Matrizen von der Berechnung der Schätzwerte abzusehen. In diesem Fall ist die Bedingung rank(L) =^! 3 nicht nur so zu verstehen, dassLvollen Rang hat, sondern auch eine genügend große Anregung der Eingänge vorliegt, um die Lösung des GleichungssystemsLθˆ = rsowie die Inverse vonLnumerisch sicher bestimmen zu können.⁹

Anstelle der numerischen Bestimmung des Ranges könnte auch überprüft werden, ob sich die Eingangs-größen genügend (linear unabhängig) geändert haben, wobei dafür geeignete Kriterien zu definieren sind.

Implementierung als RLS

Der Least-Squares-Algorithmus kann auch in einer rekursiven Form implementiert werden (siehe Ab-schnitt B.4). In dieser muss anstelle der in jedem Schritt notwendigen Invertierung der Matrix(Ψ^TΨ)⁻¹, die hier die Dimension (p+1×p+1) besitzt, bzw. der Lösung eines linearen Gleichungssystems mit p+1Unbekannten, nur eine skalare Division durchgeführt werden. Mathematisch sind die beiden Ver-fahren äquivalent, sofern die Matrizen(Ψ^TΨ)⁻¹ undPnach dem ersten Schritt identisch sind. Dies kann näherungsweise durch die Wahl sehr großer Einträge im Anfangswert P₀ oder durch die Berechnung des erstens Schritts über das „normale“ Least-Squares-Verfahren und das Fortsetzen im RLS-Verfahren erreicht werden.

9 Bei den zu dieser Untersuchung durchgeführten Simulationen und Messungen wurde die untere Grenze für die Singu-lärwerte auf1·10⁻¹⁰gesetzt, d. h. kleinere Singulärwerte wurden bei der Rangbestimmung als Null angesehen. Damit sind bei der Verwendung von Gleitkommazahlen mit 64 bit Länge keine Probleme mit der Lösung vonLθˆ=rbzw. der Inversion vonLaufgetreten.

Dem Vorteil des geringeren Rechenaufwands stehen hier aber auch Nachteile gegenüber, die durch die Notwendigkeit entstehen, Anfangswerte für θˆ und P anzugeben. Es werden dazu in [ISERMANNund MÜNCHHOF, 2011, S. 272] verschiedene Möglichkeiten angegeben. Existieren sinnvolle Schätzungen aus a-priori-Wissen (z. B. aus vorherigen Versuchen), so können diese verwendet werden.

Dies ist hier aber nicht möglich, da zum einen für θ₀ keine genaue a-priori-Schätzung angegeben wer-den kann, und zum anderen mit diesem Verfahren viele Messwerte benötigt werwer-den, um einen falschen Anfangswert zu korrigieren. Eine weitere Methode wäre,P₀=αImit großemαzu setzen. Im Grenzfall α→ ∞nähert man sich damit dem normalen LS-Verfahren an, unabhängig von dem Anfangswert θˆ₀. Aber für endlicheαwerden wieder mehrere Messwerte benötigt, um die gewünschten Werte zu erhalten.

Auch treten Probleme auf, wenn Parameter zu sekundären Eingängen erst im späteren Verlauf eines Auftrags mitgeschätzt werden können.

Letztlich bleibt nur die Möglichkeit, das RLS-Verfahren mit einem normalen LS-Verfahren zu starten.

Dies muss während eines Auftrags auch jedes Mal neu gemacht werden, wenn ein weiterer Eingang sich ändert und damit ein zusätzlicher Parameter geschätzt werden kann. Damit erhält man mit dem RLS-Verfahren immer dieselben Werte wie mit dem LS-Verfahren.

Da damit aber ohnehin Speicher für all die Daten vorgehalten werden muss, die auch für das normale LS-Verfahren benötigt werden, und auch die entsprechenden Gleichungssysteme gelößt werden müssen, wenn auch nicht so häufig, wird auf das RLS-Verfahren vollständig verzichtet.

Im Dokument Eine Methodik zur stochastischen adaptiven Qualitätsregelung (Seite 121-125)