Matroids Matheplanet Forum Index
Moderiert von matroid
Mathematik » Numerik & Optimierung » Gauss-Newton-Verfahren bei gedämpfter Schwingung
Druckversion
Druckversion
Autor
Beruf J Gauss-Newton-Verfahren bei gedämpfter Schwingung
Zaby1990
Neu Letzter Besuch: im letzten Monat
Dabei seit: 12.06.2020
Mitteilungen: 2
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Themenstart: 2020-06-16


Hallo liebes Matroids-Matheplanet-Forum,

mein derzeitiges mathematisches Problem besteht darin eine bekannte Funktion in gemessene Werte zu interpolieren, um diese analytisch beschreiben zu können. Das gemessene Signal ist eine Stromstärke in kA und entspricht einer transformierte Kondensatorentladung. Wer den genauen Hintergrund dazu lesen möchte, kann dies in der dieser Dissertation tun.

Ich habe die Formel
\(I(t)=\frac{U_0}{\omega_\mathrm{e} \cdot L}\cdot \mathrm{e}^{-\delta t}\cdot\sin{(\omega_\mathrm{e} \cdot t)}\)
mit der Eigenkreisfrequenz
\(\omega_\mathrm{e}=\sqrt{{\omega_0}^2-{\delta}^2}\)
der Resonanzkreisfrequenz
\(\omega_0 = \frac{1}{\sqrt{L\cdot C}}\)
und der Abklingkonstante
\(\delta=\frac{R}{2\cdot L}\)
gegeben. Die Formel gilt für den primären Stromkreis, meine Messdaten stammen jedoch aus dem Sekundärstromkreis. Da der Transformator nicht in Sättigung geht, kann ich die Stromstärke mit dem Übersetzungsverhältnis umrechnen. Dies wurde direkt beim Anpassen der Messdaten gemacht und muss nicht weiter beachtet werden.

Ich habe den Newton-Gauss-Algorithm nach der englischen Wikipedia in Python integriert. Zum Testen, ob die Implementierung funktioniert, habe ich zuerst das Beispiel dort abgebildet (hat funktioniert) und danach auf mein Problem erweitert.

Für die Details:

Ich habe die 3 unbekannten \(C\), \(L\) und  \(R\) im Vektor \( \vec\beta =\left(\begin{array}{c} C\\L\\R \end{array}\right) \)

Mein Residuenvektor ist \( \vec r = \left( \begin{array}{c} y_0 - I(t_0)\\y_1 - I(t_1)\\ \vdots \\y_n - I(t_n)  \end{array}    \right) \)

Die Jakobimatrix besteht aus einer Spalte je unbekannte und aus so vielen Zeilen wie Messwerten:


\( \mathbf{J} = \left(\begin{array}{rrr}
\frac{\partial r_0}{\partial C} & \frac{\partial r_0}{\partial L} & \frac{\partial r_0}{\partial R} \\
\frac{\partial r_1}{\partial C} & \frac{\partial r_1}{\partial L} & \frac{\partial r_1}{\partial R} \\
            \vdots &\vdots &\vdots \\
\frac{\partial r_n}{\partial C} & \frac{\partial r_n}{\partial L} & \frac{\partial r_n}{\partial R} \\
\end{array}\right) =
\left(\begin{array}{rrr}
\frac{\partial y_0}{\partial C} - \frac{\partial I(t_0)}{\partial C} & \frac{\partial y_0}{\partial L} - \frac{\partial I(t_0)}{\partial L} & \frac{\partial y_0}{\partial R} - \frac{\partial I(t_0)}{\partial R} \\
\frac{\partial y_1}{\partial C} - \frac{\partial I(t_1)}{\partial C} & \frac{\partial y_1}{\partial L} - \frac{\partial I(t_1)}{\partial L} & \frac{\partial y_1}{\partial R} - \frac{\partial I(t_1)}{\partial R} \\
            \vdots &\vdots &\vdots \\
\frac{\partial y_n}{\partial C} - \frac{\partial I(t_n)}{\partial C} & \frac{\partial y_n}{\partial L} - \frac{\partial I(t_n)}{\partial L} & \frac{\partial y_n}{\partial R} - \frac{\partial I(t_n)}{\partial R} \\

\end{array}\right)\)

Die Ableitungen von \(\frac{\partial y_n}{\partial C}\), \(\frac{\partial y_n}{\partial L}\), \(\frac{\partial y_n}{\partial R}\) = 0.

Die Ableitungen nach den unbekannten Parametern sind:
\(\frac{\partial I}{\partial C}(t)=-\dfrac{2\cdot U_0}{C\cdot \omega_{\mathrm{e}}\cdot\left(C\cdot R^2-4\cdot L\right)}\cdot\mathrm{e}^{-\delta\cdot t}\cdot\left(\omega_{\mathrm{e}}\cdot t \cdot\cos{\left(\omega_\mathrm{e}\cdot t\right)} - \sin{\left(\omega_{\mathrm{e}}\cdot t\right)} \right)\)

\(\frac{\partial I}{\partial L}(t)=-\dfrac{U_0}{8\cdot C\cdot L^5\cdot \omega_{\mathrm{e}}^3}\cdot\mathrm{e}^{-\delta\cdot t}\cdot\left(L\cdot\left(2\cdot L - C\cdot R^2\right)\cdot 2\cdot\omega_{\mathrm{e}} \cdot t\cdot \cos{\left(\omega_{\mathrm{e}} \cdot t\right)}+\left( 4\cdot L^2 - 4\cdot L\cdot R\cdot t + C\cdot R^3 \cdot t\right) \cdot\sin{\left(\omega_{\mathrm{e}} \cdot t\right)}\right)\)

\(\frac{\partial I}{\partial R}(t)=-\frac{U_0}{8\cdot C \cdot L^4 \cdot  \omega_\mathrm{e}^3}\cdot\mathrm{e}^{-\delta t}\cdot \left( C \cdot L \cdot R \cdot  2\cdot \omega_\mathrm{e} \cdot \cos{\left( \omega_\mathrm{e} \cdot t \right)}
+ \left( 4\cdot L\cdot t -2\cdot C \cdot R\cdot L - C \cdot R^2 \cdot t  \right) \cdot \sin{\left( \omega_\mathrm{e} \cdot t \right)} \right)\)

Diese Ableitungen habe ich in die Jakobimatrix eingesetzt. An das zusätzliche Vorzeichen (wegen \(0 - \frac{\partial I(t_n)}{\partial L}\) etc.) habe ich gedacht.

In jedem Berechnungsschritt werden zunächst  \(\vec r\) und \(\mathbf{J}\) berechnet. Damit werden die unbekannten neu bestimmt durch
\(\vec\beta_\mathrm{neu} = \vec\beta + (\mathbf{J}^\intercal \cdot \mathbf{J})^{-1}\cdot \mathbf{J}^\intercal \cdot \vec r \),

Aus anderen Ermittlungsansätzen (s.o. genannte Dissertation) kenne ich recht genaue Werte für die Maschinenparameter
\(C= 0.008 F\)
\(L= 0.00031 H\)
\(R= 0.162 R\)
die mir als Schätzwert der Parameter für den ersten Iterationsschritt dienen.

Aufgrund der Maschinencharakteristik entspricht die gemessene Stromkurve jedoch nur bis kurz nach dem Strommaximum der oben angegebenen Formel. Deshalb reduziere ich die Daten bis zu dem Punkt, wo die theoretische mit der gemssenen Kurve noch gut übereinstimmt.

Die Implementierung habe ich wie folgt gelöst (zum Testen kann ich gern die CSV-Datei irgendwie zuschicken):
python
import matplotlib.pyplot as plt
import numpy as np
import glob
import pandas as pd
from matplotlib.cm import get_cmap
 
fig, (ax1, ax2) = plt.subplots(2)
df = pd.read_csv('M5_122.csv')
print(df['U_C'].iloc[0])
 
 
c_real = 0.008
l_real = 0.000315
r_real = 0.182
 
 
def omega0(c, l):
    return (1 / np.sqrt(l * c))
 
 
def delta(l, r):
    return (r / (2 * l))
 
 
def omega_e(omega0, delta):
    return np.sqrt(omega0**2 - delta**2)
 
 
def I_ke(t, u0, c, l, r):
    _omega0 = omega0(c, l)
    _delta = delta(l, r)
    _omega_e = omega_e(_omega0, _delta)
 
    ret = (u0 / (_omega_e * l)) * \
        (np.exp(-1 * _delta * t)) * np.sin(_omega_e * t)
    return ret
 
 
def residual(t, I, u0, c, l, r):
    return I - I_ke(t, u0, c, l, r)
 
 
def dIdL(t, u0, c, l, r):
    _omega0 = omega0(c, l)
    _delta = delta(l, r)
    _omega_e = omega_e(_omega0, _delta)
 
    term1 = (-1) * (u0 / (8 * (c) * (l**5) * (_omega_e**3)))
 
    term2 = np.exp(-1 * _delta * t)
 
    term3 = (l * (2 * l - c * r**2) * 2 * _omega_e * t * np.cos(_omega_e * t) +
             (4 * l**2 - 4 * l * r * t + c * r**3 * t) * np.sin(_omega_e * t))
    return term1 * term2 * term3
 
 
def dIdC(t, u0, c, l, r):
    _omega0 = omega0(c, l)
    _delta = delta(l, r)
    _omega_e = omega_e(_omega0, _delta)
 
    term1 = (-2) * (u0 / (_omega_e * c * (c * r**2 - 4 * l)))
 
    term2 = np.exp(-1 * _delta * t)
 
    term3 = _omega_e * t * np.cos(_omega_e * t) - np.sin(_omega_e * t)
 
    return term1 * term2 * term3
 
 
def dIdR(t, u0, c, l, r):
    _omega0 = omega0(c, l)
    _delta = delta(l, r)
    _omega_e = omega_e(_omega0, _delta)
 
    term1 = (-1) * (u0 / (8 * (c) * (l**4) * (_omega_e**3)))
 
    term2 = np.exp(-1 * _delta * t)
 
    term3 = (c * l * r * 2 * _omega_e * np.cos(_omega_e * t) +
             (4 * l * t - 2 * c * r * l - c * r**2 * t) * np.sin(_omega_e * t))
 
    return term1 * term2 * term3
 
# Berechnung mit den bekannten Daten
df['I_calc'] = I_ke(df['t'], df['U_C'].iloc[0], c_real, l_real, r_real)
 
# Begrenzen der Messdaten
df.drop(range(df['I'].argmax()+6,df['t'].argmax()+1),inplace=True)
 
Anzahl_unbekannter = 3
 
R = np.ones((len(df['t']), 1))
Jacobi = np.ones((len(df['t']), Anzahl_unbekannter))
 
c_guess = 0.008  # in F
l_guess = 0.000315  # in H
r_guess = 0.122  # in Ohm
 
beta = np.array([[c_guess], [l_guess], [r_guess]])
 
diff = 100
k = 0
sum_residuals_old = 0
sumres = []
while (diff > 0.00001) & (k < 10**4):
    k += 1
    print(k)
    for i in range(len(df['t'])):
        R[i, 0] = residual(df['t'].iloc[i], df['I'].iloc[i], df['U_C'].iloc[0], beta[
                           0], beta[1], beta[2])
 
    sum_residuals = np.sum(np.square(R))
    diff = np.abs(sum_residuals - sum_residuals_old)
    sum_residuals_old = sum_residuals
    sumres.append(sum_residuals)
    for i in range(np.shape(Jacobi)[0]):
        Jacobi[i, 0] = (-1) * dIdC(df['t'].iloc[i],
                                   df['U_C'].iloc[0], beta[0], beta[1], beta[2])
        Jacobi[i, 1] = (-1) * dIdL(df['t'].iloc[i],
                                   df['U_C'].iloc[0], beta[0], beta[1], beta[2])
        Jacobi[i, 2] = (-1) * dIdR(df['t'].iloc[i],
                                   df['U_C'].iloc[0], beta[0], beta[1], beta[2])
 
    beta -= np.dot(np.dot(np.linalg.inv(
        np.dot(Jacobi.transpose(), Jacobi)), Jacobi.transpose()), R)
 
    print(beta)
 


Zur Verdeutlichung gibt es auch noch ein Bild, damit ihr euch besser vorstellen könnt, wovon ich eigentlich rede.





Nun zum Problem:
Nach wenigen Durchläufen werden alle Parameter zu NANs. Und ich weiß nicht, woran es liegt. Ist das Verfahren ungeeignet? Wenn ja, gibt es ein Verfahren, dass dafür geeignet ist? Leider dachte ich immer, dass ich die Mathematik eh nicht mehr brauche und war deshalb in den Vorlesungen wenig aufmerksam.

Fällt jemanden vielleicht ein Fehler auf in den Ableitungen oder sonstigem? (Ich weiß, das zu prüfen ist sehr aufwendig, aber vielleicht hat jemand Freude daran.)

Ich hoffe, ich habe euch genügend Informationen gegeben ohne euch zu erschlagen. Ich freue mich auf eure Antworten und Ideen.

Viele Grüße und vielen Dank für eurer Interesse
Zaby



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
DetlefA
Aktiv Letzter Besuch: im letzten Monat
Dabei seit: 12.11.2007
Mitteilungen: 83
Aus: Berlin-Tegel
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.1, eingetragen 2020-07-13


Hallo Zaby,

ich habe verstanden, dass Du die Parameter R,L,C einer gedämpften Sinusschwingung aus Messwerten bestimmen möchtest indem Du das nichtlineare Gleichungssystem mit Hilfe der Jacobimatrix löst.

Zur Bestimmung der Parameter einer gedämpften Schwingung ist es nicht nötig ein nichtlineares Gleichungssystem zu lösen, das läßt sich mit linearer Algabra direkt erledigen.

Die Vorgehensweise habe ich hier

dargelegt und hier

implementiert.
Ich hoffe, das ist nützlich
Cheers
Detlef



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Zaby1990
Neu Letzter Besuch: im letzten Monat
Dabei seit: 12.06.2020
Mitteilungen: 2
Zum letzten BeitragZum nächsten BeitragZum vorigen BeitragZum erstem Beitrag  Beitrag No.2, vom Themenstarter, eingetragen 2020-07-13


Hallo Detlef,

vielen Dank für deine Antwort. Ich habe mir das angesehen und werde das auf jeden Fall testen und mit meinem Ergebnis vergleichen. Ich sehe allerdings das Problem, dass ich Eigenkreisfrequenzen <600Hz habe, die Messungen aber mit 200kHz oder sogar 1MHz Abtastrate durchgeführt werden.


Ich habe diese Woche meinen Fehler gefunden: es haben sich Fehler in den Ableitungen eingeschlichen. Mit den korrekten Ableitungen komme ich mit 5-10 Itertationsschritten (abhängig von den Schätzwerten) auf sehr gute Lösungen. Ist die Schätzung zu schlecht, konvergiert das Problem nicht.

Die richitgen Ableitungen sind:

\[\frac{\partial I}{\partial C}(t)= -\dfrac{U_0}{2 \cdot \omega_{e}^3 \cdot L^2}\cdot \mathrm{e}^{-\delta\cdot t}\cdot
\left[ \frac{1}{C^2} \cdot \Bigl(\omega_e \cdot t\cdot \cos{\left(\omega_e\cdot t\right)}-\sin{\left(\omega_e\cdot t\right)}\Bigr) \right]\] \[\frac{\partial I}{\partial L}(t)= -\dfrac{U_0}{2 \cdot \omega_{e}^3 \cdot L^2}\cdot \mathrm{e}^{-\delta\cdot t}\cdot \left[\left( \frac{1}{C\cdot L}-\frac{R^2}{2\cdot L^2}\right) \cdot \Bigl(\omega_e \cdot t\cdot \cos{\left(\omega_e\cdot t\right)}-\sin{\left(\omega_e\cdot t\right)}\Bigr) + 2\cdot \omega_{e}^2 \cdot (1-\delta\cdot t)\cdot \sin{\left(\omega_e\cdot t\right)}       \right]\] \[\frac{\partial I}{\partial C}(t)= -\dfrac{U_0}{2 \cdot \omega_{e}^3 \cdot L^2}\cdot \mathrm{e}^{-\delta\cdot t}\cdot
\left[ \frac{1}{C^2} \cdot \Bigl(\omega_e \cdot t\cdot \cos{\left(\omega_e\cdot t\right)}-\sin{\left(\omega_e\cdot t\right)}\Bigr) \right]\]


Nochmals vielen Dank für die Hilfe, auch von denen, die mir nur privat geantwortet haben.

Viele Grüße
Zaby



Eine Notiz zu diese Forumbeitrag schreiben Notiz   Profil  Quote  Link auf diesen Beitrag Link
Zaby1990 hat die Antworten auf ihre/seine Frage gesehen.
Zaby1990 hat selbst das Ok-Häkchen gesetzt.
Zaby1990 wird per Mail über neue Antworten informiert.
Neues Thema [Neues Thema]  Druckversion [Druckversion]

 


Wechsel in ein anderes Forum:
 Suchen    
 
All logos and trademarks in this site are property of their respective owner. The comments are property of their posters, all the rest © 2001-2020 by Matroids Matheplanet
This web site was originally made with PHP-Nuke, a former web portal system written in PHP that seems no longer to be maintained nor supported. PHP-Nuke is Free Software released under the GNU/GPL license.
Ich distanziere mich von rechtswidrigen oder anstößigen Inhalten, die sich trotz aufmerksamer Prüfung hinter hier verwendeten Links verbergen mögen.
Lesen Sie die Nutzungsbedingungen, die Distanzierung, die Datenschutzerklärung und das Impressum.
[Seitenanfang]