Parameterschätzung, überbestimmtes LGS, Kalman-Filter

[Randy B. (rbrecker)]

Hallo zusammen,

ich habe folgendes Problem:

Einem Experiment liegt der folgende, simple lineare Zusammenhang 
zugrunde:

m1(t) = a * m2(t) + b * m3(t)

Hierbei sind m1, m2, m3 Messwerte, a und b die Kenngrößen des Systems, 
also konstant und zu ermitteln.

Mit

m1(t1) = a * m2(t1) + b * m3(t1)
m1(t2) = a * m2(t2) + b * m3(t2)

habe ich ein LGS, was sich wegen der störungsbehafteten Messwerte m1, 
m2, m3 nur näherungsweise lösen lässt.

m1, m2, m3 lassen sich kontinuierlich über der Zeit messen, ich bekomme 
also eine beliebig lange Folge von Messwerten. Und schritthaltend mit 
den Messwerten möchte ich nun die Parameter a und b schätzen.

Wonach ich nun suche, ist ein iterativer Algorithmus, bei dem ich 
jeweils durch Eingabe von [m1, m2, m3](t) eine neue Schätzung für [a, 
b](t) bestimme (vllt. im Sinne von LMS). Also so ähnlich wie ein 
exponentieller Mittelwertschätzer.

Die unterschiedlichen Streuungen der Messungen kann man auch zunächst 
außer Acht lassen.

Irgendwie "riecht" das nach Kalman-Filter, jedoch stecke ich da 
irgendwie fest.

[Ob S. (Firma: 1984now) (observer)]

Randy B. schrieb:

> ich habe folgendes Problem:
[...]

> Wonach ich nun suche, ist ein iterativer Algorithmus, bei dem ich
> jeweils durch Eingabe von [m1, m2, m3](t) eine neue Schätzung für [a,
> b](t) bestimme (vllt. im Sinne von LMS). Also so ähnlich wie ein
> exponentieller Mittelwertschätzer.

Grundsätzlich: Fortgesetzte lineare Regression.

Und dann ein wenig mit den Ergebnissen spielen. Eine genaue Vorschrift 
kann man nicht angeben, weil die von der zunkünftigen Entwicklung 
abhängt, die man natürlich zu keinem Zeitpunkt t kennen kann (jedenfalls 
nicht in unserem Universum).

Alles, was man tun kann, ist: aus der Historie einen Trend ableiten und 
beten, dass das System den ermittelten Trend zumindest in der näheren 
Zukunft in etwa beibehält. Kann es tun, kann es aber auch nicht tun...

[Joe G. (feinmechaniker)]

Wie verhalten sich deine Messwerte m1, m2 und m3?
Wenn diese jeweils nur (normalverteilt) um ihren Mittelwert schwanken, 
wird es schwierig. Dann hast du nämlich nur EINE Gleichung mit 2 
Unbekannten.

[Randy B. (rbrecker)]

Joe G. schrieb:
> Wie verhalten sich deine Messwerte m1, m2 und m3?
> Wenn diese jeweils nur (normalverteilt) um ihren Mittelwert schwanken,
> wird es schwierig. Dann hast du nämlich nur EINE Gleichung mit 2
> Unbekannten.

Danke für den Hinweis! Ich bin wohl auf dem Holzweg ...

Es handelt sich um eine Messung an einem DC-Motor:

m1(t) ~ Ue(t)
m2(t) ~ i(t)
m3(t) ~ rpm(t)

bestimmt werden sollen a ~ Ri und b ~ Km des Motors.
Und solange ich die Belastung des Motors nicht ändern kann, werden die 
einzelnen Messung wohl (fast) linear abhängig sein, so dass man keine 
brauchbaren Ergebnisse aus der Schätzung bekommt.

Ich stelle die Frage ggf. nochmal in einem anderen Thread, wie man mit 
einem PWM-Steller, bei dem man den Motorstrom in der PWM(ON)-Phase genau 
messen kann, ggf. doch auf a ~ Ri schließen kann. Und dies mit einem 
Testlauf des Motors ohne Last.

V/ESC (für BLDC-Motoren) kann ja auch alle notwendigen Motorparameter 
messen ...

[Joe G. (feinmechaniker)]

Du könntest ja bei gleicher Last auf zwei unterschiedlichen 
Arbeitspunkten Messen, also z.B. zwei unterschiedliche Ue(t) bei 
gleicher Last. Somit werden es zwei Gleichungen mit zwei Unbekannten und 
die können gelöst werden (siehe Anhang)

[Johann (Gast)]

Randy B. schrieb:
> Es handelt sich um eine Messung an einem DC-Motor:
>
> m1(t) ~ Ue(t)
> m2(t) ~ i(t)
> m3(t) ~ rpm(t)

Nur um es richtig zu verstehen: Wir reden hier von einem stationären 
Betrieb, jedoch alle Größen sind hier abhängig von der Zeit. Falls 
dynamisch, müsste man auch die Schwungmasse berücksichtigen

Ui = Uklemme - Ri*I    (=innere Spannung)
Ui = kv * Omega
Mi = km * I     (=inneres Moment)
dOmega/dt * Inertia = Mi - Mlast

Für sehr schnelle Änderungen sollte man dann noch die Ankerinduktivität 
berücksichtigen.

[Randy B. (rbrecker)]

Johann schrieb:
> Nur um es richtig zu verstehen: Wir reden hier von einem stationären
> Betrieb, jedoch alle Größen sind hier abhängig von der Zeit. Falls
> dynamisch, müsste man auch die Schwungmasse berücksichtigen
>

> Ui = Uklemme - Ri*I    (=innere Spannung)
> Ui = kv * Omega
> Mi = km * I     (=inneres Moment)
> dOmega/dt * Inertia = Mi - Mlast
>

>
> Für sehr schnelle Änderungen sollte man dann noch die Ankerinduktivität
> berücksichtigen.

Ja

[A. S. (rava)]

du möchtest also das hier:

https://de.wikipedia.org/wiki/Multiple_lineare_Regression#Sch%C3%A4tzung_des_Parametervektors_mit_der_Kleinste-Quadrate-Sch%C3%A4tzung

aber in iterativ?

Spontan fällt mir dazu aufgrund der Matrixinversen keine effiziente 
Umformung ein.

Wenn ich jetzt schnell etwas zusammenhacken müsste, würde ich versuchen 
ein Fenster von Messwerten zu behalten und damit in jedem Tick die 
Parameter neu zu schätzen - eine art "windowed LMSE". Wenn diese Werte 
dann noch zu sehr rauschen, kann man entweder das Fenster zu vergrößern 
oder die Schätzwerte in einem postprocessing-Schritt zu glätten. Dann 
folgt das Ergebnis natürlich nichtstationärem Verhalten langsamer.

[Detlef _. (detlef_a)]

Hi,

Du hast ein überbestimmtes Gleichungssystem m1=a*m2+b*m3
Mit den m1(1)...m1(n),m2(1)...m2(n),m3(1)...m3(n) als Spaltenvektoren 
bekommst Du das überbestimmte lineare Gleichungssystem

M * [a]= [m2(1) m3(1)]*[a]=[m1(1)]
    [b]  [m2(2) m3(2)] [b] [m1(2)]
          .......            ....

Die lineare Regression als Lösung dazu lautet
[a]     = inv(MT * M)  MT  m1
[b]

mit inv als Matrixinversion und MT als Transponierter von M.
MT*M ist eine 2x2Matrix mit den Elementen
Summealler m2*m2
Summealler m3*m3
und Summealler m2*m3

MT*m1 ist ein 2er Spaltenvektor mit den Elementen
Summealler m2*m1 und Summealler m3*m1

Schätzwerte für diese 5 Summen berechnest Du jetzt mit einem 
exponentiellen Mittelwertschätzer.

Ein exponentieller Mittelwertschätzer geht so:
x(k+1)=x(k)+Eingangswert-(1/c)*x(k)
Wenn x sich nicht mehr ändert gilt x(k+1)=x(k), also 
Eingangswert=(1/c)*x(k)

Der Schätzer liefert ein gemitteltes Signal der Eingangswerte, das um 
den Faktor c größer ist. c bestimmt wie 'schnell' das Ding ist.

Und für die Summen m2*m2, m3*m3, m2*m3, m2*m1 und m3*m1 nimmst Du dann 
jeweils so einen exponentiellen Schätzer. Das c bestimmt, wie schnell 
die Lösung für a,b den Messwerten folgt.

Das kostet wenig Aufwand, insbesondere wenn c eine Potenz von 2 ist und 
geht sehr gut.

been there, done that
math rulez!
Cheers
Detlef

[Detlef _. (detlef_a)]

Hier nochmal mit lesbarer Formatierung

Die Formel:

$$\begin{align*} m_1(k)=a \cdot m_2(k)+b \cdot m_3(k) \end{align*}$$

Als Matrizengleichung

$$\begin{bmatrix} m_2(1) & m_3(1) \\ m_2(2) & m_3(2) \\ \dots \end{bmatrix} \cdot \begin{bmatrix} a \\ b \\ \end{bmatrix} = \begin{bmatrix} m_1(1) \\ m_1(2) \\ \dots \end{bmatrix}$$

$$\begin{bmatrix} m{_2}^T & \ m{_3}^T \\ \end{bmatrix} \cdot \begin{bmatrix} a \\ b \\ \end{bmatrix} = \begin{bmatrix} m{_1}^T \\ \end{bmatrix} = M \cdot \begin{bmatrix} a \\ b \\ \end{bmatrix}$$

Die Lösung:

$$\begin{bmatrix} a \\ b \\ \end{bmatrix} = (M^T \cdot M)^{-1} \cdot M^T \cdot m_1^T \newline = \begin{bmatrix} \sum_k m_2 \cdot m_2 & \sum_k m_2 \cdot m_3 \\ \sum_k m_2 \cdot m_3 & \sum_k m_3 \cdot m_3 \\ \end{bmatrix}^{-1} \cdot \begin{bmatrix} \sum_k m_2 \cdot m_1 \\ \sum_k m_3 \cdot m_1 \\ \end{bmatrix}$$

Die fünf verschiedenen Summen der Produkte werden mit exponentiellen 
Mittelwertbildung kontinuierlich mitgezogen

$$x_{k+1}=x_k+Eingangswert- (1/c) \cdot x_k$$

Für

$$x_{k+1}=x_k$$

ändert sich nix mehr und es gilt

$$Eingangswert = (1/c) \cdot x_k$$

Das Ding liefert also einen Schätzwert, der Faktor c zu gross ist.
Je größer c, umso langsamer ist der plot.

math rulez.
Cheers
Detlef

[Randy B. (rbrecker)]

Detlef _. schrieb:
> Die fünf verschiedenen Summen der Produkte werden mit exponentiellen
> Mittelwertbildung kontinuierlich mitgezogen
> xk+1=xk+Eingangswert−(1/c)⋅xk
> x_{k+1}=x_k+Eingangswert- (1/c) \cdot x_k

Danke für die Ausführungen, das ist eigentlich genau das, was ich 
gesucht habe!

Den exp. Mittelwert kenne ich so:

$$x_{k+1} = c \cdot x_k + (1 - c) \cdot v$$

mit v als neuem Wert.

[Christoph M. (mchris)]

>exp. Mittelwert kenne ich so:

Für mich neu, dass man das exponentieller Mittelwert nennt. Für mich 
beschreibt das das Verhalten eines RC-Tiefpasses bezw. eines PT1-Gliedes 
in der Regelungstechnik.

[Joe G. (feinmechaniker)]

Leider läßt sich das Problem nicht so lösen, wie es von Detlef 
beschrieben wurde.
Es bleibt ja immer noch eine Gleichung mit zwei Unbekannten. Wenn das 
Problem, wie von Detlef skizziert, als lineares Gleichungssystem in 
Matrixform formuliert wird, stößt man sehr schnell auf ein Problem mit 
der inversen Matrix inv(MT * M).
Letztlich muß dabei die Determinante von MT*M gebildet werden. Bei 
kleinen Messfehlern der Messvektoren m2 und m3 geht diese gegen null. 
Bei idealen Messwerten von m2 und m3 wird sie genau null.
Das bedeutet praktisch: Je geringer die Messfehler, um so größer der 
Schätzfehler von a und b – sehr ungünstig! Im Anhang mal die komplette 
Rechnung, wenn es interessiert.

[Detlef _. (detlef_a)]

Hi,

die Determinante ist die Energie von m2 mal der Energie von m3 minus 
Quadrat der Kreukorrelation zwischen m2 und m3. Die Determinante ist nur 
Null für m2=m3.

Habs mal getestet, geht wunderbar.

math rulez!
Cheers
Detlef

clear
n=10000;
m2=randn(1,n);
m3=randn(1,n);
a=1.12345;
b=2.2345;
m1=a*m2+b*m3;

M=[m2.' m3.'];
M.'*M
[m2*m2.' m2*m3.' ; m2*m3.' m3*m3.']
coff=inv(M.'*M)*M.'*m1.'

m22=zeros(1,n);
m33=zeros(1,n);
m23=zeros(1,n);
m21=zeros(1,n);
m31=zeros(1,n);
ar=zeros(1,n);
br=zeros(1,n);
c=256;
for(k=2:n)
  m22(k)=m22(k-1)+m2(k)*m2(k)-(1/c)*m22(k-1);
  m33(k)=m33(k-1)+m3(k)*m3(k)-(1/c)*m33(k-1);
  m23(k)=m23(k-1)+m2(k)*m3(k)-(1/c)*m23(k-1);
  m21(k)=m21(k-1)+m2(k)*m1(k)-(1/c)*m21(k-1);
  m31(k)=m31(k-1)+m3(k)*m1(k)-(1/c)*m31(k-1);
  M=[m22(k) m23(k) ; m23(k) m33(k)]*(1/c);
  coff=inv(M.'*M)*[m21(k) ; m31(k)]*1/c;
  ar(k)=coff(1);
  br(k)=coff(2);
end;

return

[Joe G. (feinmechaniker)]

Hallo Detlef,
mit viel Rauschen geht es, jetzt nehme mal reale Werte eines 
Gleichstrommotors mit relativ wenig Rauschen, oder lasse mal das 
Rauschen ganz weg (konstante Werte) dann wird die Matrix auch in deiner 
Rechnung sigulär.

clear
n=100;

m2=ones(1, n) * 0.15 + 0.001*randn(1,n);
m3=ones(1, n) * 930 + 2.3*rand(1,n);

a=0.4;
b=0.006387;

m1=a*m2+b*m3;

M=[m2.' m3.'];
M.'*M
[m2*m2.' m2*m3.' ; m2*m3.' m3*m3.']
coff=inv(M.'*M)*M.'*m1.'

m22=zeros(1,n);
m33=zeros(1,n);
m23=zeros(1,n);
m21=zeros(1,n);
m31=zeros(1,n);
ar=zeros(1,n);
br=zeros(1,n);
c=256;
for(k=2:n)
  m22(k)=m22(k-1)+m2(k)*m2(k)-(1/c)*m22(k-1);
  m33(k)=m33(k-1)+m3(k)*m3(k)-(1/c)*m33(k-1);
  m23(k)=m23(k-1)+m2(k)*m3(k)-(1/c)*m23(k-1);
  m21(k)=m21(k-1)+m2(k)*m1(k)-(1/c)*m21(k-1);
  m31(k)=m31(k-1)+m3(k)*m1(k)-(1/c)*m31(k-1);
  M=[m22(k) m23(k) ; m23(k) m33(k)]*(1/c);
  coff=inv(M.'*M)*[m21(k) ; m31(k)]*1/c;
  ar(k)=coff(1);
  br(k)=coff(2);
end;

Anmerkung zu den Parametern:
m2 = 150 mA
m3 = 930 1/s
m1 = 6V

a = 0.4 Ohm
b = 0.006387 Wb

[Detlef _. (detlef_a)]

Ja klar. Wenn Du immer das gleiche mißt kannst Du keine 2 Konstanten 
bestimmen.

[Matthias M. (matthias_m222)]

Beim RLS-Algorithmus braucht man zB keine Matrix-Inversion. Das Beispiel 
dort -> https://de.m.wikipedia.org/wiki/RLS-Algorithmus erklärt es ganz 
gut.

[Joe G. (feinmechaniker)]

Vielen Dank für den Hinweis zum RLS-Algorithmus. Für den Parameter b 
(Motorkonstante) liefert er sehr gute Ergebnisse, der Parameter a 
(Motorwiderstand) gibt ein brauchbares Ergebnis.

@author: Joe G.
"""
import numpy as np
N = 10000
a = 0
b = 0

# RLS Algorithm parameters
lambda_ = 1  # Forgetting factor
delta = 1e-8    # Small positive constant to initialize P

# Initialization
n_params = 2  # Number of parameters (a and b)
P = np.eye(n_params) / delta
theta = np.zeros(n_params)  # Initial guess for [a, b]

# Simulated data (example data)
m1 = np.random.normal(6, 0.05, N)  # Motor voltage
m2 = np.random.normal(0.15, 0.01, N)  # Motor current
m3 = np.random.normal(930, 5, N)  # Angular speed

# Stack m2 and m3 as feature matrix X
X = np.vstack((m2, m3)).T

# RLS algorithm loop
for i in range(len(m1)):
    x = X[i, :].reshape(-1, 1)  # Current feature vector
    y = m1[i]  # Current target value
    
    # Prediction error
    y_pred = np.dot(theta, x)
    e = y - y_pred
    
    # Kalman gain vector
    k = P @ x / (lambda_ + x.T @ P @ x)
    
    # Parameter update
    theta = theta + k.flatten() * e
    
    # Covariance matrix update
    P = (P - k @ x.T @ P) / lambda_

# Output the estimated parameters
a, b = theta
print("Estimated a:", a)
print("Estimated b:", b)