Parameterschätzung, überbestimmtes LGS, Kalman-Filter

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 eine neue Schätzung für a, b 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...

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.

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 ...

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 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

1

> Ui = Uklemme - Ri*I    (=innere Spannung)

2

> Ui = kv * Omega

3

> Mi = km * I     (=inneres Moment)

4

> dOmega/dt * Inertia = Mi - Mlast

5

>

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

Ja

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.

Hi,

Du hast ein überbestimmtes Gleichungssystem m1=am2+bm3 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. MTM ist eine 2x2Matrix mit den Elementen Summealler m2m2 Summealler m3m3 und Summealler m2m3

MTm1 ist ein 2er Spaltenvektor mit den Elementen Summealler m2m1 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 m2m2, m3m3, m2m3, m2m1 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

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

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.

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.

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.

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=am2+bm3;

M=[m2.' m3.']; M.'M [m2m2.' m2m3.' ; m2m3.' 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

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.001randn(1,n); m3=ones(1, n) * 930 + 2.3rand(1,n);

a=0.4; b=0.006387;

m1=am2+bm3;

M=[m2.' m3.']; M.'M [m2m2.' m2m3.' ; m2m3.' 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

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

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.

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.

1

@author: Joe G.

2

"""

3

import numpy as np

4

N = 10000

5

a = 0

6

b = 0

7

8

# RLS Algorithm parameters

9

lambda_ = 1  # Forgetting factor

10

delta = 1e-8    # Small positive constant to initialize P

11

12

# Initialization

13

n_params = 2  # Number of parameters (a and b)

14

P = np.eye(n_params) / delta

15

theta = np.zeros(n_params)  # Initial guess for [a, b]

16

17

# Simulated data (example data)

18

m1 = np.random.normal(6, 0.05, N)  # Motor voltage

19

m2 = np.random.normal(0.15, 0.01, N)  # Motor current

20

m3 = np.random.normal(930, 5, N)  # Angular speed

21

22

# Stack m2 and m3 as feature matrix X

23

X = np.vstack((m2, m3)).T

24

25

# RLS algorithm loop

26

for i in range(len(m1)):

27

    x = X[i, :].reshape(-1, 1)  # Current feature vector

28

    y = m1[i]  # Current target value

29

30

    # Prediction error

31

    y_pred = np.dot(theta, x)

32

    e = y - y_pred

33

34

    # Kalman gain vector

35

    k = P @ x / (lambda_ + x.T @ P @ x)

36

37

    # Parameter update

38

    theta = theta + k.flatten() * e

39

40

    # Covariance matrix update

41

    P = (P - k @ x.T @ P) / lambda_

42

43

# Output the estimated parameters

44

a, b = theta

45

print("Estimated a:", a)

46

print("Estimated b:", b)