Kennlinien-Linearisierung mit pow()

[bounceboy (Gast)]

Guten Tag,

ich versuche mich gerade daran, eine stromgesteuerte Bremse mit einem 
AVR zu anzusteuern (funktioniert mit PWM auch ganz gut). Der AVR soll 
zudem über ein Display das aufgebrachte Drehoment anzeigen. Da der 
Drehmomentverlauf nicht ganz linear ansteigt, habe ich das Drehmoment an 
der Welle über den Duty Cycle meiner PWM gemessen und den Verlauf in 
Excel mit einer Trendlinie nachvollzogen. Laut Excel lautet die Formel 
zur Trendlinie 6. Grades:

y = 1E-10x^6 - 3E-08x^5 + 4E-06x^4 - 0,0003x^3 + 0,0087x^2 - 0,0227x

Mithilfe dieser Formel soll der AVR nun das Drehmoment berechnen:

unsigned int pwm_prozent;               // PWM-Prozentwert 0...100
double moment_nm;

double berechne_moment_nm(void)
{
moment_nm = (((double)(0.0000000001 * (pow(pwm_prozent, 6)))) - 
((double)(0.00000003 * (pow(pwm_prozent, 5)))) + ((double)(0.000004 * 
(pow(pwm_prozent, 4)))) - ((double)(0.0003 * (pow(pwm_prozent, 3)))) + 
((double)(0.0087 * (pow(pwm_prozent, 2)))) - ((double)(0.0227 * 
pwm_prozent)));

return moment_nm;

}

Allerdings liefert mir der Controller nicht die erwarteten Ergebnisse.
Als Controller verwende ich einen Mega16 mit dem Pollin-Board, als 
Compiler dient CodeVision, die math.h habe ich ebenfalls eingebunden!

Kann das mit meinen Datentypen zusammenhängen? Habe auch schon mal 
gehört, dass es Probleme mit der pow-Funktion gibt, kann mir hier jemand 
weiter helfen?

Vielen Dank schon mal im Voraus!

[Karl H. (kbuchegg)]

Schon mal was von der Horner Regel gehört?

y = a*x^2 + b*x

ist identisch zu

y = ( ( a * x ) + b ) * x;

(x wird einfach herausgehoben).
Damit wirst du die ganzen pow Aufrufe los, was der Genauigkeit
deiner Rechnerei zugute kommen dürfte (und schneller geht es
nebenbei auch noch)

y = (((((( 1E-10 * x ) - 3E-08 ) *  x + 4E-06 ) * x - 0,0003 ) * x +
      0,0087 ) * x - 0,0227 ) * x;

Wie sinnvoll es ist, die 6., 5. und 4. Potenz mit einzubeziehen
ist eine andere Frage, die man untersuchen müsste. Mit einem
WinAVR double, hast du gerade mal 6 bis 7 signifikante Stellen.
Und das wird dann schon mächtig eng.

[yalu (Gast)]

> Wie sinnvoll es ist, die 6., 5. und 4. Potenz mit einzubeziehen ist
> eine andere Frage,

Der Wertebereich für x geht von 0 bis 100. Für x=100 liegen alle fünf
Summanden des Terms in einer ähnlichen Größenordnung, so dass man
keinen weglassen sollte.

[yalu (Gast)]

Die Funktion sieht für eine Bremse ja schon recht seltsam aus :)
Oder sollte x vielleicht doch nur von 0 bis 1 gehen?

[I_ H. (i_h)]

pow scheint im AVR-GCC nicht zu funktionieren. Für ganzzahlige 
Exponenten ist das aber schnell selber implementiert.

[Jörg W. (dl8dtl) (Moderator)]

I_ H. wrote:

> pow scheint im AVR-GCC nicht zu funktionieren.

Was genau funktioniert nicht?

[I_ H. (i_h)]

Garnix funktioniert. Das Problem hatte ich schonmal vor einer Weile. 
Spuckt irgendwelche Ergebnisse aus, die mit der Realität nix zu tun 
haben.

Beitrag "pow()-Funktion"

[Patrick N. (pneuberger)]

Hallo,

bau die Formel um. Nimm das Horner-Schema von Karl-Heinz. Weiterhin 
würde ich darüber nachdenken, die doubles in der Formel mit 1E10 zu 
multiplizieren, mit unsigned Ganzzahlen zu rechnen und das Ganze am 
Schluß zu dividieren, weil double beim GCC nur float-Genauigkeit hat.


MfG

Patrick

[Patrick N. (pneuberger)]

Waaah... signed meine ich natürlich.

[Jörg W. (dl8dtl) (Moderator)]

I_ H. wrote:
> Garnix funktioniert. Das Problem hatte ich schonmal vor einer Weile.
> Spuckt irgendwelche Ergebnisse aus, die mit der Realität nix zu tun
> haben.
>
> Beitrag "pow()-Funktion"

Schnee von gestern, davon abgesehen, dass dort auf Threads von 2004
verwiesen wird.  Kurzes Testprogramm:

Enter x, y: 3 3
pow(3, 3) = 27
Enter x, y: 3 5
pow(3, 5) = 243
Enter x, y: 5 3
pow(5, 3) = 125
Enter x, y: 2.54 3
pow(2.54, 3) = 16.3871
Enter x, y: 1e-5 2
pow(1e-05, 2) = 1e-10
Enter x, y: 1e-5 25
pow(1e-05, 25) = 0
Enter x, y: 2.5 2.5
pow(2.5, 2.5) = 9.88211
Enter x, y: -2.5 3
pow(-2.5, 3) = -15.625
Enter x, y: 10 3
pow(10, 3) = 1000
Enter x, y: 3 10
pow(3, 10) = 59049
Enter x, y: 2 10
pow(2, 10) = 1024
Enter x, y: 9 9
pow(9, 9) = 3.8742e+08
Enter x, y:

Ziemlich viel für "garnix funktioniert", oder?

Welche Zahlen soll ich dir noch eintippen?

[I_ H. (i_h)]

Also vor ein paar Monaten hatte ich immernoch das Problem, auch mit 
einer aktuellen Verion. Btw. eintippen - du hast das aber schon auf 'nem 
AVR laufen lassen?

[Jörg W. (dl8dtl) (Moderator)]

I_ H. wrote:

> Btw. eintippen - du hast das aber schon auf 'nem
> AVR laufen lassen?

Ja, na klar.  Ich habe immer irgendwo ein Gerippe für eine einfache
stdio-Applikation via UART rumliegen, da habe ich einfach ein
scanf() und ein printf() in einer Endlosschleife reingestopft.

Wie gesagt, wenn du für die aktuelle Version (1.6.x) irgendetwas
findest, wo pow() nicht geht (oder nicht im Rahmen der machbaren
Genauigkeit), dann schreibe einen Bugreport.  Momentan sehe ich
keinen offenen.

[I_ H. (i_h)]

Hab damals das neueste draufgemacht was der gentoo portage tree hergab 
(öhm... libc 1.4.6), und damit hat's nicht funktioniert. Hatte damals 
auch den passenden Bugreport gefunden, und mir dann pow(10, ...) selber 
geschrieben.

Aber gut, wenn's inzwischen gefixt ist. Inzwischen ist 1.6.x auch im 
portage tree angekommen.

[Jörg W. (dl8dtl) (Moderator)]

Ja, seit dem Weggang von Henrik Brix Andersen ist Gentoo's Support
für die AVR-Toolchain ziemlich hinten runtergefallen, habe ich den
Eindruck.

[andreas (Gast)]

Hallo zusammen!

An bounceboy: Bist du sicher, dass deine Interpolationsfunktion stimmt. 
Ich hatte mal den Fall bei einem NTC. Habe Werte aufgenommen und wollte 
die von Excel interpoliert haben. Ab dem dritten Grad kam nur noch 
Schrott raus. Die Interpolationspolynome haben die Stützstellen nicht 
mehr getroffen!

Gib doch einmal die Interpolationsfunktion in Excel ein, und lass sie 
über die Messwerte legen.

[bounceboy (Gast)]

Hallo, vielen Dank erst mal für eure Antworten...des Rätsels Lösung lag 
doch etwas näher als gedacht. Die Funktion, die mir Excel ausgegeben hat 
ist völliger Schwachsinn, warum auch immer?!? Habe das ganze mal in 
meinen grafikfähigen Taschenrechner eingedaddelt und komme auf das 
gleiche wie Yalu (Danke an dieser Stelle an Yalu und Andreas). Gibt es 
irgendwelche andere Möglichkeiten, Kennlinien mit durch Punkte 
mathematisch anzunähern?

[yalu (Gast)]

Eine Interpolation bzw. Approximation durch ein Polynom n-ten Grades
ist nicht immer geschickt. Oft ist es so, dass ein niedriger Grad zu
ungenau wird, da die Kurve zu weit an den vorgegebenen Punkte
vorbeiläuft, und ein hoher Grad zwar die Punkte trifft, aber zu
unschönen Wellen zwischen den Punkten führt.

Am besten ist es, wenn man die Physik hinter der ganzen Sache
verstanden hat und eine allgemeine Funktion aufstellen kann, deren
Parameter dann mit Hilfe der Stützpunkte geschätzt werden können.

Oder man erkennt durch scharfes Anschauen der Punkte und ein wenig
Ausprobieren, welcher Funktionstyp gut passen würde und schätzt die
Parameter dann ebenfalls aus den Stützpunkten.

An deiner Stelle würde ich die Messpunkte einfach mal posten,
vielleicht hat ja jemand eine passende Idee.

[bounceboy (Gast)]

Dann stell ich mal hier meine Parameter ein:

PWM[%]   M[Nm]
  0      0
  5      0,10
 10      0,36
 15      0,86
 20      1,32
 25      1,89
 30      2,29
 35      2,74
 40      2,96
 45      3,34
 50      3,59
 55      3,92
 60      4,19
 65      4,36
 70      4,57
 75      4,67
 80      4,86
 85      4,99
 90      5,04
 95      5,14
100      5,19

Wäre nett, wenn mir hier jemand mit einer sinnvollen und hinreichend 
genauen Funktion für den Wertebereich von 0 bis 100% aushelfen könnte.
Vielen Dank

[Εrnst B. (ernst)]

Nimm die Stützstellen, wie in deinem Post angegeben, und interpolier 
linear zwischen denen. Sollte hinreichend genau sein, und ist mit einer 
Zeile C erledigt.

[yalu (Gast)]

Ich habe gerade mal ein wenig rumgespielt. Mit Polynomen bin ich nicht
sehr weit gekommen, mit

$$y = a\left(1-\frac 1 {1+bx^2}\right)$$

bzw.

y = a - a / (1 + b * x*x)

schon eher. Durch "visuelles Probieren" bin ich auf a=6 und b=6.5e-4
gekommen.

Die Funktion ist relativ leicht zu berechnen (auf jeden Fall leichter
als ein Polynom 6ten Grades), und der Approximationsfehler könnte noch
halbwegs in der Größenordnung der Messfehler liegen.

Die Ergebnisse im Vergleich (ym: gemessen, ya: approximiert):

  x   ym    ya   Fehler
  0  0.00  0.00   0.000
  5  0.10  0.10  -0.004
 10  0.36  0.37   0.006
 15  0.86  0.77  -0.094
 20  1.32  1.24  -0.082
 25  1.89  1.73  -0.157
 30  2.29  2.21  -0.075
 35  2.74  2.66  -0.080
 40  2.96  3.06   0.099
 45  3.34  3.41   0.070
 50  3.59  3.71   0.124
 55  3.92  3.98   0.057
 60  4.19  4.20   0.014
 65  4.36  4.40   0.038
 70  4.57  4.57  -0.004
 75  4.67  4.71   0.041
 80  4.86  4.84  -0.023
 85  4.99  4.95  -0.043
 90  5.04  5.04   0.002
 95  5.14  5.13  -0.014
100  5.19  5.20   0.010

[Bernd (Gast)]

Kleiner Tip zu Excel:

y = 1E-10x^6 - 3E-08x^5 + 4E-06x^4 - 0,0003x^3 + 0,0087x^2 - 0,0227x

Das Ergebnis 1E-10x^6 ist zu ungenau. Dazu die errechnet Funktion mit 
der rechten Maustaste anklicken und als Zahlenformat Wissenschaftlich 
mit 3 oder 4 Nachkommastellen wählen.

Benötigst Du nicht die Umkehrfunktion?

Du möchtest z.B. mit 3 Nm bremsen und benötigst die Pulsweite dazu.
X und Y vertauschen und dann nochmal eine neue Funktion suchen.

Oder y = a - a / (1 + b * x*x) nach x umstellen

x = sqr((a/(a-y)-1)/b)

[Bernd (Gast)]

Das Polynom laut Excel funktioniert ganz gut:

y = 5,437E-03x6 + 1,890E-01x5 - 2,976E+00x4
  + 1,503E+01x3 - 3,162E+01x2 + 3,667E+01x

Aber
y = sqr((a/(a-x)-1)/b)
mit a = 5.989, b = 0.00065 kostet weniger Rechenleistung.

Im Wertebereich x >= 0 und x <= 5.19 ist die Funktion stetig.
Aber x muss vor dem Aufruf der Wurzelfunktion überprüft werden.

int nm2pwm(float nm)
{
int pwm;

  if (nm < 0)
    pwm = 0;
  else
  if (nm > 5.19)
    pwm = 100;
  else
    pwm = sqr((5.989/(5.989-nm)-1)/0.00065); 

  return pwm;
}

Eine weitere Möglichkeit wäre, die Funktion aus 3 einfachen Polynomen 2. 
Ordung in der Form: y = a0 + a1*x + a2*x*x zusammenzusetzten.
Die 3 Bereiche wären dann z.B. 0 - 1 Nm, 1 - 4 Nm und 4 - 5.19 Nm

int nm2pwm(float nm)
{
  if (nm < 0)
    return 0;
  else
  if (nm < 1)
    return  // hier Polynom 1
  else
  if (nm < 4)
    return  // hier Polynom 2
  else
  if (nm < 5.19)
    return  // hier Polynom 3
  else
    return 100;
}

Etwas Mühe macht es, die Polynome and den Grenzstellen genau 
aufeinandertreffen zu lassen. Aber das Ergebnis wäre schnell und genau.