Forum: Digitale Signalverarbeitung / DSP Approximation fuer PT1


Announcement: there is an English version of this forum on EmbDev.net. Posts you create there will be displayed on Mikrocontroller.net and EmbDev.net.
von Dergute W. (derguteweka)


Bewertung
0 lesenswert
nicht lesenswert
Moin,

Ich hab' ein paar gemessenen Temperaturen eines konstant beheizten Oinks 
ueber ein paar Zeiten:
z.B. sowas:
          0       :  23
          2       :  28
          4          30
          7          34
         10          37
         15          40
         20          44
         30          48
         45          53
         50          54
         60          55
         75          57
         90          58
        105          58


Hat jemand grad zufaellig einen wenigzeiler fuer GNU Octave parat, wie 
ich aus diesen Daten mir die Groessen A und Z rausapproximieren kann, 
wenn ich glaub', dass das Oink ein PT1 Verhalten hat - also sowas von 
dem Kaliber:

Temp(zeit)=Temp(0)+A*(1-exp(-zeit/Z)

Gruss
WK

von chris (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Mach doch den quadratischen Fehler und variiere A und Z zufällig und 
nimm den Treffer mit dem kleinsten Fehler. Ansonsten kannst Du es auch 
mit Gradientenabstieg versuchen

von Josef (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Mit Octave kann ich nicht dienen, aber hier eine Loesung
mit maxima.
load(lsquares);

M : matrix( [0,23],[2,28],[4,30],[7,34],[10,37],[15,40],[20,44],[30,48],[45,53],[60,55],[75,57],[90,58],[105,58] );

lsquares_estimates( M, [x,y], y = B + A * (1-exp(-x/Z)), [B,A,Z]);


ergibt:

[[B = 24.50205474148759, A = 33.71562608547086, Z = 23.47286576998589]]

von Josef (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Natuerlich mit Startwert temp(0)=23:

lsquares_estimates( M, [x,y], y = 23 + A * (1-exp(-x/Z)), [A,Z]);

ergibt:

[[A = 34.73653966932802, Z = 21.17559226183299]]

von Dergute W. (derguteweka)


Bewertung
0 lesenswert
nicht lesenswert
Moin,

Josef schrieb:
> [[A = 34.73653966932802, Z = 21.17559226183299]]

Ha, sogar schon vorgekaut, danke!

Gruss
Wolfram

von Josef (Gast)


Bewertung
1 lesenswert
nicht lesenswert
Bei Octave laeuft das wohl unter lsqcurvefit.

Ein Beispiel mit einer Exponentialfunktion ist in der Doku von Octave.

von Detlef _. (detlef_a)


Bewertung
0 lesenswert
nicht lesenswert
Da brauchts' keine schwere Artillerie wie von Josef gezeigt, das kann 
man on the fly mit jedem Mikrocontroller und linearer Algebra 
erschlagen.

>>>>>Temp(zeit)=Temp(0)+A*(1-exp(-zeit/Z)

A krieg ich raus 34.59 und Z 23.76.

math rulez!
Cheers
Detlef


clear
w=[ 0 23 2 28 4 30 7  34 10 37 15 40 20 44 30 48 .....
45 53 50 54 60 55 75  57 90 58 105 58 ];

t=w(1:2:end);
temp=w(2:2:end);
plot(t,temp,'r.-')

endwert=58.5;
y=(endwert-23)-(temp-23);
y=y(:);
x=t(:);

x=x(1:end);
y=y(1:end);

n=length(y);
M=[x(:) ones(n,1)];
coff=inv(M'*M)*M'*log(y);
a=coff(2);
b=coff(1);
plot(x,log(y),x,M*coff);
xr=x;
yr=exp(a)*exp(b.*xr);
yr=23-(yr-(58.1-23));
semilogy(x,y,'b.-',xr,yr,'r.-')

plot(t,temp,xr,yr)
return

von Sigi (Gast)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Detlef _. schrieb:
> Da brauchts' keine schwere Artillerie wie von Josef gezeigt, das kann
> man on the fly mit jedem Mikrocontroller und linearer Algebra
> erschlagen.
>
>>>>>>Temp(zeit)=Temp(0)+A*(1-exp(-zeit/Z)
>
> A krieg ich raus 34.59 und Z 23.76.

Nein, um Gottes Willen! Das Problem ist (leider, wie so oft)
nicht linear, d.h. nicht linear in den zu optimierenden
Parametern A und Z. Wenn's z.B. um die Form A*exp(-z+Z) gehen
würde, dann könnte es einfach auf eine lineare Form reduziert
werden (und z.B. per PseudoInverse invert(M^T*M)*M^T das
Optimum bestimmen). Man hätte aber immer noch grössere
numerische Probleme.

Vlt. eine kleine Analogie: du willst ein Integral per Aussummieren
grob berechnen, änderst dann aber die Skala von x nach exp(x).
Dadurch ändern sich ja die "Boxen" unter der Kurve, und zwar
in X-Richtung. Die Summe entspricht damit überhaupt nicht mehr
dem Integral. Beim obigen Problem ist es genauso. Durch deine
"Linearisierung" änderst du die Topologie im L2 extrem stark.
Sie ist nicht mehr isotrop und nicht mehr ortsunabhängig.

Ich habe mal ein einfaches Programm geschrieben, dass A und Z
variiert und je Paarung den Fehler ausrechnet. Die Bilder im
Anhang zeigen dann die Fehler, zur besseren Visualisierung
habe ich den Logarithmus des Fehlers geplottet. Bild 1 zeigt
den groben Überblick, Bild2 einen feineren Ausschnitt als
Kontourplot. Das so bestimmte Optimum ist A=34.68 und Z=21.236,
deckt sich also mit der Lösung aus Octave.

Man könnte per Verfeinerung das Optimum noch besser eingrenzen,
aber in Anbetracht der schlechten Eingangsbedingung (Addition
mit Temp(0) !!, ich würd's umformulieren) spielt das kaum eine
Rolle.

von Detlef _. (detlef_a)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
>>>>>
Nein, um Gottes Willen! Das Problem ist (leider, wie so oft)
nicht linear, d.h. nicht linear in den zu optimierenden
Parametern A und Z.
<<<<<<<<

Das ist nicht das Optimum, ist mir bekannt. Seit 200Jahren lösen 
Ingenieure das Fitting einer Exponentialfunktion mit der Logarithmierung 
und damit Reduktion auf ein lineares Problem.

Aktuell machen das auch die Seuchenforscher genauso, die logarithmieren 
und legen ne Gerade in die Fallzahlen.

Das hat noch nie nicht funktioniert.

Ein linearer Fit geht auf einem uC, eine nichtlineare Optimierung nicht.

Im Übrigen ist meine Optimierung, obwohl linear, bisschen besser als 
Deine. Das dürfte eigentlich nicht sein.

math rulez
Cheers
Detlef


clear
w=[ 0 23 2 28 4 30 7  34 10 37 15 40 20 44 30 48 .....
45 53 50 54 60 55 75  57 90 58 105 58 ];

t=w(1:2:end);
temp=w(2:2:end);
plot(t,temp,'r.-')

endwert=58.5;
y=(endwert-23)-(temp-23);
y=y(:);
x=t(:);

x=x(1:end);
y=y(1:end);

n=length(y);
M=[x(:) ones(n,1)];
coff=inv(M'*M)*M'*log(y);
a=coff(2);
b=coff(1);
plot(x,log(y),x,M*coff);
xr=x;
yr=exp(a)*exp(b.*xr);
yr=23-(yr-(58.1-23));
semilogy(x,y,'b.-',xr,yr,'r.-')


zr=23+34.68*(1-exp(-x/21.236))

plot(t,temp,'b.-',xr,yr,'r.-',xr,zr,'k.-')
grid
title(sprintf('Sigis Abweichungen %.2f meine Abweichungen %.2f ',....
    sum(abs(temp-zr.')),sum(abs(temp-yr.'))))
xlabel('Zeit [sec]')
ylabel('Temperatur [°C]')
legend('Messwerte', 'Ick','Sigi')
return

von Sigi (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Zu meinem "Optimum" muss natürlich auch gesagt werden,
dass es kein Optimum im math. Sinne ist. Ich suche ja
nur das Optimum in einem 2D-Array mit endlicher
Auflösung.

Zu deiner Lösung: ich habe mir einen Maier gesucht, bis
ich auf deine 9.9 statt 14.78 als Fehlerwert kam. Für
meine Berechnungen hatte ich nur die TO-Formel von ganz
Oben verwendet, du verwendest dagegen

yr=exp(a)*exp(b.*xr);
yr=23-(yr-(58.1-23));

Ausserdem nimmst du als Formel für den Fehler die L1-Norm,
das verändert alles (ich habe nach L2 optimiert). Wenn ich
mein Array-Ansatz (was ja mit Mathe nichts zu tun hat, ist
nur eine Holzhammermethode) gegen unterschiedliche Normen
anwende (L1, L2 und L-Unendlich), dann komme ich auf
folgende Tabelle mit Fehlerwerten:

-------------------------------------------------------
TO-Formel: z = temp[0] + A*(1-exp(-x/Z))
-------------------------------------------------------
User:       L1              L2              Linf
-------------------------------------------------------
IKE :       14.782785       22.886193       2.2075594
IKE2:       9.9540454       11.449178       1.6975611
-------------------------------------------------------
L1  :       9.6679087       12.787407       2.0401627
L2  :       10.379489       10.966322       1.8822150
Linf:       11.959534       13.987099       1.6848860
-------------------------------------------------------

IKE: dein Optimum mit TO-Formel
IKE2: dein Optimum mit deiner Formel
I1: mein Optimum bzgl. L1-Norm
I2: mein Optimum bzgl. L1-Norm
Iinf: mein Optimum bzgl. L-Unendlich-Norm

Bei meinen 3 Optima sieht man dann sehr schön,
das je nach Norm auch ein "gutes" Optimum gefunden
wird, dessen Fehler auch besser ist als die anderen.

Das dein Fehler bei L1 trotzdem noch besser ist
als meiner liegt an deiner Formel, die von der
TO-Formel abweicht. Prüf mal bitte dein
Optimum-Paar a,z mit der TO-Formel, dann müsstest
du wie ich auf 14.78 kommen.

von Sigi (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Detlef _. schrieb:
> Ein linearer Fit geht auf einem uC, eine nichtlineare Optimierung nicht.

Mit meiner Array-Methode wäre ein uC schnell
überfordert, ich würde aber obiges Optimum
z.B. per Gradientenmethode bestimmen. Und
das geht, je nach Formulierung, sogar
schneller als über den linearen Ansatz.

von Detlef _. (detlef_a)


Bewertung
0 lesenswert
nicht lesenswert
Hi Sigi,

ja, alles mit der heissen Nadel gestrickt bei mir, das war ein hack.

Zwei Sachen möchte ich sagen:

Erstens: Exponentialfkt. mit Logarithmierung und Fit mit Pseudoinversen 
funktioniert immer, auf dem kleinsten uC.

Zweitens: Iterationen nach Qualitätskriterien ( Gradientenverfahren, 
Echocanceller, Channel Equalizer) zicken immer, dauert 8 Wochen bis der 
Sch läuft.

Been there, done that. War so, ist so und wird solange bleiben, bis die 
Mathematiker mal bessere nichtlineare Optimierungsverfahren basteln. 
Kann dauern ;)).

math rulez, anyhow.
Cheers
Detlef

von Sigi (Gast)


Angehängte Dateien:

Bewertung
0 lesenswert
nicht lesenswert
Detlef _. schrieb:
> Zweitens: Iterationen nach Qualitätskriterien ( Gradientenverfahren,
> Echocanceller, Channel Equalizer) zicken immer, dauert 8 Wochen bis der
> Sch läuft.

Was'n Blödsinn: Erstens, ich habe gerade 2 Verfahren
implementiert. Beides PlainVanilla Gradientenabstieg,
einmal nach L2 und einmal nach L_unendlich. Das
Runterschreiben hat jeweils 10-20 Minute gedauert,
dazu nochmal je 10-20 Minuten Testen.
Zweitens, schau dir die Zielfunktion an (T0-A*(1-...)),
die ist visuell extremst gutmütig. Ich habe nur die Ableitungen
ausgerechnet, von den aktuellen Vektor (a,z) abgezogen
und die Trajektorie beobachtet. Wie erwartet (war mir
99.9% sicher), kein gezicke, schnelle Konvergenz (in
Bezug auf Rechenzeit als auch auf Anzahl Iterationsschritte).
Ob das auf einem uC läuft oder nicht ist total schnuppe,
war auch nicht gefordert, Ziel ist nur das Optimum.
Du kannst jetzt auch noch ergänzende Verfahren wie
adjungierte Gradienten oder 2teOrdungsverfahren
(konvergiert noch schneller als das Adjungierte), und
du kannst dich dabei im ggs zu den Ingenieursbasteleien
darauf verlassen, dass sie auch von Profinumerikern
verifiziert, publiziert, ebenfalls von Profis anerkannt
und weltweit millionenfach getestet wurden, auch wieder
von Profis.

Detlef _. schrieb:
> Erstens: Exponentialfkt. mit Logarithmierung und Fit mit Pseudoinversen
> funktioniert immer, auf dem kleinsten uC.

Du wirst zwar noch immer glauben, ich sei nicht dazu
in der Lage, das ebenfalls so zu machen (in der Medizin
etc. kann man so für eine schnelle Visualisierung vorgehen,
für numerisch anspruchsvolle Geschichten aber nicht zu
gebrauchen). Ich habe mal im Anhang eine Graphik hochgeladen,
die deutlichst zeigt, warum. Zu sehen ist eine EXP-Kurve
mit zwei Fehlerkurven (je +0.01 bzw -0.01 hinzuaddiert)
und davon dann den Logarithmus geplottet. Du siehst ja
hoffentlich ganz schön, dass der Fehler nach Hinten raus
mit O(exp(x)) bzw. O(-exp(x)) wächst. In z.B. der Medizin
kann das ja toleriert werden, sei's zur Visualisierung
oder nur zur groben Abschäzungen für was weiss ich.
Und damit nun zur schlechten Konditionierung des Problems,
speziell bezogen auf die Addition mit Temp(0): Hier wackelt
nicht der Schwanz mit dem Hund, sondern die Schwanzspitze.
Der Fehler in Temp(0) überträgt sich auch das ganze Intervall.
Und mit den Fehlern +/-O(exp(x)) bei deinem Verfahren hast
du dir einen zweiten Schwanz an den Hund gebastelt, der auch
noch stark am Wackeln ist. Schon das erste Problem würde
mich sofort veranlassen, die Approximationsfunktion geeigneter
zu wählen, aber ist ja nur meine bescheidene Meinung.
Deine Implementierung hat noch weitere Probleme bzw. Fehler,
die suchst du nun aber bitte selbst.

Detlef _. schrieb:
> Been there, done that. War so, ist so und wird solange bleiben, bis die
> Mathematiker mal bessere nichtlineare Optimierungsverfahren basteln.
> Kann dauern ;)).

Also du optimierst (mit allen bereits diskutierten und
wirklich leicht nachzuvollziehenden Fehlern) im L2,
bemasst dein Ergebnis dann im L1 (weiss du überhaupt,
was hier L1 macht, und wenn schon anderes Mass, dann
BITTE L_unendlich, macht im ggS zu L1 wenigstens Sinn,
bei all der Sinnlosigkeit in L(k1) zu optimieren und
in L(k2) zu testen, k1!=k2) und sprichst dann bei
Mathematikern von basteln. Ich bin selbst kein
Mathematiker, habe mich nur in verschiedenen Bereichen
vertieft, uA auch in Numerik (I+II+nur0.25*III
wegen uniteresanter Spezalisierung), da wird nichts
gebastelt, also spar dir dein Fachbashing.
Wenn du ein besseres Verfahren haben willst:
1. statt invert(M^T*M)*M^T verwende doch mal
invert(M^T*V*M)*M^T*W mit W*W=W^T*W=V und V als
Diagonalmatrix. Diese einfache Form ergibt sich, wenn
in L2 nicht das cannonische Innenprodukt, sondern
das gewichtete Innenprodukt verwendet wird. Die
Diagonalelemente von V sind dann die Gewichte,
die Diagonalelemente von W sind gleich den Wurzeln
der Diagonalelemente von V. Diese Gewichte kannst
du dann so wählen, dass sie den zweiten von mir
beschriebenen Fehler ausgleichen. Leite das doch einfach
mal aus der Optimierungsformel per Ableitung nach
dem Lösungsvektor ab, dann kriegst du auch ein bischen
Übung, kann nie schaden, braucht man immer.
2. Rechne damit dann eine Lösung aus, berechne die
Differenz zum Ziel und Approximiere dann mit der selben
Formel die Differenz. Das Ergebnis wird dann zum ersten
Ergebnis hinzuaddiert und liefert dann wegen der
Additivität ein wesentlich besseres Ergebnis. Diesen
Schritt kannst du auch noch ein bis zweimal wiederholen,
womit du dein ein sauschnelles und hochgenaues
Iterationsverfahren hast (dazu gibt's viele Paper).

Kann dauern: Du kennst Bücher wie z.B. Numerical Recipies?
Oder überhaupt ein Buch über Numerik (habe hier den
Hammerlin, für praktische Sachen besser als der Schwarz),
dazu am Besten noch ein Funktionalanalysisbuch (habe
hier die Einführungsbibel von Heuser, damit sind dann
alle Fragen zu L(k) und allen damit verbundenen Ungleichungen
erschlagen, braucht man wirklich immer)). Wenn du damit
nicht fündig wirst..

Und da Mathematik auch immer mit Geschichte passenderweise
zusammengelehrt wird:
Detlef _. schrieb:
> Seit 200Jahren lösen
> Ingenieure das Fitting einer Exponentialfunktion mit der Logarithmierung
> und damit Reduktion auf ein lineares Problem.

So aus dem Kopf, Gauss hatte vor ca 200 seine Methode
der kleinsten Quadrate entwickelt. Im ggs zu Cauchy
war aber Gauss bzgl. Publizieren bekanntich eine faule
Socke. Es hätte also so schon lange gedauert, bis die
Methode gross Einzug in andere Wissenschafstbereiche
gefunden hätte. Dazu kommt, dass diese Methoden erst
durch die moderne Numerik grössere Verbreitung fand,
und die ist grossteils erst in der zweiten Hälfte des
letzten Jahrhunderts entstanden (ich weiss, Numerik
gab's auch schon früher).

von Sigi (Gast)


Bewertung
0 lesenswert
nicht lesenswert
Und weil ich gerade die Maxima-Seite auf habe:
Das Ergebnis meiner Optimierung ist ein wenig
(eiglich deutlich) besser als das von Maxima.
Das mag eine von wenigen Ausnahmen sein, aber
hier sieht man dann doch ein, warum Profitools
bzw. Bibliotheken für z.B. Mathematika etc.
so gottverdammt teuer sind, eben wegen ihrer
Zuverlässigkeit. Einige dieser Tools bzw. Libs
kosten schnell 4-5stellig im Jahr an Lizenz.
Den Ärger, den man sich damit aber einspart
kompensiert das aber total.

Antwort schreiben

Die Angabe einer E-Mail-Adresse ist freiwillig. Wenn Sie automatisch per E-Mail über Antworten auf Ihren Beitrag informiert werden möchten, melden Sie sich bitte an.

Wichtige Regeln - erst lesen, dann posten!

  • Groß- und Kleinschreibung verwenden
  • Längeren Sourcecode nicht im Text einfügen, sondern als Dateianhang

Formatierung (mehr Informationen...)

  • [c]C-Code[/c]
  • [avrasm]AVR-Assembler-Code[/avrasm]
  • [code]Code in anderen Sprachen, ASCII-Zeichnungen[/code]
  • [math]Formel in LaTeX-Syntax[/math]
  • [[Titel]] - Link zu Artikel
  • Verweis auf anderen Beitrag einfügen: Rechtsklick auf Beitragstitel,
    "Adresse kopieren", und in den Text einfügen




Bild automatisch verkleinern, falls nötig
Bitte das JPG-Format nur für Fotos und Scans verwenden!
Zeichnungen und Screenshots im PNG- oder
GIF-Format hochladen. Siehe Bildformate.

Mit dem Abschicken bestätigst du, die Nutzungsbedingungen anzuerkennen.