Nukix
Messages postés5Date d'inscriptionjeudi 2 juillet 2009StatutMembreDernière intervention15 juillet 2009
-
7 juil. 2009 à 14:51
alusipdranreb
Messages postés2Date d'inscriptionjeudi 22 avril 2010StatutMembreDernière intervention25 mai 2010
-
25 mai 2010 à 09:09
Bonjour à tous !
je suis nouveau sur le fofo, on m'a dit qu'il y avait des pros qui trainaient dans le coin, alors je suis venu voir.
Voilà,
je suis ingénieur et j'essaye d'utiliser un programme de regression polynomiale pour trouver un polynome approchant la position d'un ballon stratosphèrique en fonction du temps car la précision du capteur gps ne me convient pas.
Voici ce que j'aimerai approximer:
il se trouve que le programme me sort une approximation mais pas de tous les points, il s'arrête avant la fin.
Qu'en pensez vou ?
j'avoue que j'ai du mal à comprendre le programme.
(Le code est en bas du message)
Public x() As Double 'Données d'entrée (xi et yi)
Public y() As Double
Public ai_p() As Double 'Valeurs des coefficients du polynôme de tendance
Sub main()
Dim i As Integer, j As Integer
Dim ligne As String
Dim virg_deb As Long
Dim virg_fin As Long
Dim a As String
i = 0
'Boucle de récupération des données xi et yi à traiter
Open "valeurs_init.txt" For Input As #1
While Not EOF(1)
i = i + 1
ReDim Preserve x(1 To i)
ReDim Preserve y(1 To i)
Line Input #1, ligne
'Remplacement des blancs par des séparateurs (chr(47) = "/")
ligne = Replace(ligne, Chr(9), Chr(47))
ligne = Replace(ligne, Chr(32), Chr(47))
ligne = Chr(47) & ligne & Chr(47)
'Récupération des données x(i) et y(i)
'Recherche des séparateurs (/) et stockage des valeurs indiquées entre les séparateurs
virg_deb = InStr(1, ligne, Chr(47))
virg_fin = InStr((virg_deb + 1), ligne, Chr(47))
x(i) = CDbl(Mid(ligne, (virg_deb + 1), (virg_fin - virg_deb - 1)))
'Appel de la procédure de régression polynômiale :
'le premier argument est le degré du polynôme de tendance chosi, x() et y() sont les données à traiter
'ai_p() est un argument de retour renvoyant les coef ai du polynôme
Call regression_polynomiale(12, x(), y(), ai_p())
Dim xf() As Double
Dim yf() As Double
ReDim xf(1 To UBound(x))
ReDim yf(1 To UBound(y))
'Boucle générant un fichier texte contenant les valeurs xf(i) et yf(i) approchées
Open "valeurs_fin.txt" For Output As #1
For i = 1 To UBound(x)
xf(i) = x(i)
yf(i) = ai_p(0)
For j = 1 To (UBound(ai_p()))
yf(i) = yf(i) + ai_p(j) * xf(i) ^ j
Next j
Print #1, xf(i) & Chr(9) & Format(yf(i), "####0.000")
Next i
End Sub
Sub regression_polynomiale(deg_pol As Integer, xi() As Double, yi() As Double, coef_pol() As Double)
Dim i As Integer, j As Integer
Dim k As Integer, l As Integer
Dim nb_pts As Integer
Dim mat_A() As Double, mat_B() As Double
Dim mat_At() As Double, mat_Bt() As Double
Dim dim_mat As Integer
dim_mat = deg_pol + 1
'Nombre de points
nb_pts = UBound(xi)
If nb_pts <> UBound(yi) Then MsgBox "Problème, il n'y a pas le même nombre de valeurs pour les x et les y"
ReDim mat_A(1 To dim_mat, 1 To dim_mat)
ReDim mat_B(1 To dim_mat)
ReDim mat_At(1 To dim_mat, 1 To dim_mat)
ReDim mat_Bt(1 To dim_mat)
ReDim coef_pol(0 To deg_pol)
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'' Etape 1 :
'' Mise en équation du problème de minimisation par la méthode des moindres carrés
'' on cherche à minimiser : E(a0,..., ap) = somme(yi-somme(aj*xi^j)²) où p est le degré du pôlynome
'' les dérivées partielles de E par rapport à chacun des coef du polynôme de tendance doivent être nulles,
'' c'est une condition nécessaire !
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' les xi et yi sont transmis en paramètres
' la procédure renvoie alors en sortie 2 matrices mat_A et mat_B tq :
' mat_A * vect_des_aj = mat_B
Call mise_en_éq(nb_pts, deg_pol, xi(), yi(), mat_A(), mat_B())
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'' Etape 2 :
'' On réalise alors la triangulation de la mat_A et les opérations équivalentes sur mat_B
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
' mat_A et mat_B sont données en arguments d'entrée
' mat_At et mat_Bt sont renvoyées en arguments de sortie, après triangulation de mat_A
Call triangulation(dim_mat, mat_A(), mat_B(), mat_At(), mat_Bt())
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'' Etape 3 :
'' Cette procédure résout alors le système et renvoie les coefficients aj du pôlynome de tendance
''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''''
'les matrices mat_At et mat_Bt sont données en arguments d'entrée
'les coef aj sont renvoyées en arguments de sortie sous le nom de variable coef_pol
Call resolution(dim_mat, mat_At(), mat_Bt(), coef_pol())
End Sub
Sub mise_en_éq(n As Integer, p As Integer, vect1() As Double, vect2() As Double, mat1() As Double, mat2() As Double)
Dim i As Integer, j As Integer, k As Integer
For k = 0 To p ' on balaie les lignes
For j = 0 To p ' on balaie les colonnes
mat1(k + 1, j + 1) = 0
For i = 1 To n ' on somme tous les xi^(k+j)
mat1(k + 1, j + 1) = mat1(k + 1, j + 1) + vect1(i) ^ (k + j)
Next i
Next j
Next k
For k = 0 To p ' on balaie les lignes
mat2(k + 1) = 0
For i = 1 To n ' on somme tous les (yi*xi^(k))
mat2(k + 1) = mat2(k + 1) + (vect2(i) * vect1(i) ^ (k))
Next i
Next k
End Sub
Sub triangulation(n As Integer, mat1() As Double, mat2() As Double, mat3() As Double, mat4() As Double)
Dim i As Integer, j As Integer, k As Integer
Dim mat3k() As Double, mat4k() As Double
ReDim mat3k(1 To n, 1 To n, 1 To n) As Double
ReDim mat4k(1 To n, 1 To n) As Double
For i = 1 To n
mat4k(i, 1) = mat2(i)
For j = 1 To n
mat3k(i, j, 1) = mat1(i, j)
Next j
Next i
For k = 1 To (n - 1)
For i = 1 To n
For j = 1 To n
If i <= k Then
mat3k(i, j, (k + 1)) = mat3k(i, j, k) ' on conserve les mêmes valeurs
mat4k(i, (k + 1)) = mat4k(i, k)
Else
mat4k(i, (k + 1)) = mat4k(i, k) - (mat3k(i, k, k) * mat4k(k, k) / mat3k(k, k, k))
If j <= k Then
mat3k(i, j, (k + 1)) = 0
Else
mat3k(i, j, (k + 1)) = mat3k(i, j, k) - (mat3k(i, k, k) * mat3k(k, j, k) / mat3k(k, k, k))
End If
End If
Next j
Next i
Next k
For i = 1 To n
mat4(i) = mat4k(i, n)
For j = 1 To n
mat3(i, j) = mat3k(i, j, n)
Next j
Next i
End Sub
Sub resolution(n As Integer, mat3() As Double, mat4() As Double, coef_pol() As Double)
Dim i As Integer, j As Integer, k As Integer
Dim xsol() As Double
Dim som_aijxj As Double
ReDim xsol(1 To n)
For i = (n - 1) To 1 Step (-1)
som_aijxj = 0
For j = i + 1 To n
som_aijxj = som_aijxj + mat3(i, j) * xsol(j)
Next j
xsol(i) = 1 / mat3(i, i) * (mat4(i) - som_aijxj)
coef_pol(i - 1) = xsol(i)
Next i
marinmarais
Messages postés104Date d'inscriptionlundi 11 avril 2005StatutMembreDernière intervention16 juillet 20101 10 juil. 2009 à 09:19
Bonjour,
Si je comprends bien ton problème, tu as un ballon sonde doté d'un récepteur GPS. A l'intant t, tu enregistres la date et l'altitude de ton ballon sonde.
Et c'est là que je ne suis pas sur de comprendre. Tu cherches à interpoler la position de ton ballon n'importe quand entre ton instant zéro et ta dernière mesure, où alors tu cherches à savoir où est ton ballon entre ta dernière mesure et la suivante que tu ne connais pas encore.
Pour la première solution, je ne m'embêterais même pas avec les moindres carrés pour chercher un polynôme. Je ferai une FFT. On s'en fout des fréquences et des longueurs d'onde. Tu auras simplement une fonction numérique correpsondant à tous tes points de mesures. Quelle que soit l'allure de ta courbe en fonction de la météo, etc, ton calcul fonctionnera, alors que tel ou tels polynômes avec les moindres carrés...
Pour la deuxième solution, à partir de tes deux dernières mesures, je calculerais la dernière tendance de ta courbe. Une approximation linéaire en attendant le dernier enregistrement.
marinmarais
Messages postés104Date d'inscriptionlundi 11 avril 2005StatutMembreDernière intervention16 juillet 20101 10 juil. 2009 à 12:15
OK Florent, j'ai pigé.
Dans ce cas là, avant tout, je me demande si tu ne devrais pas utiliser un autre logiciel de calcul. Tu sais utiliser MatLab ? Tu pourrais résoudre ton problème en 5 min avec ce soft si tu es relativement familier avec des calculs matriciel usuels.
Je vois toujours deux solutions :
- La FFT. Tu appliques une FFT avec toutes tes mesures. Si tu veux limiter lisser ton signal, et bien tu ne prends pas en compte les plus petites fréquences calculées. Si tu dois programmer ça en VB, tu vas passer des plombes. Commande MatLab : Y=FFT(X)... Ton résultat est dans l'ensemble des complexes mais tu as peu de boulot pour obtenir le résultat réel de la forme :
altitude(t)= a0 + [somme de 1 à n](a_i*sin(omega_i*t+Phase_i)
Tu pourras aisément dériver ton résultat pour étudier les forces en jeu et déduire ton coefficient de frottement.
- La méthode des moindres carrés : Tu ne connais pas, a priori, le degré optimal de ton polynôme. Je te conseille, avant d'implémenter une solution générale pour tous tes ballons, d'effectuer plusieurs ajustements en étudiant les résidus et surtout, en regardant à partir de quel degré de polynôme ceux-ci suivent la loi normale. C'est la méthode que j'ai adoptée pour optimiser le choix de polynômes de calibration de capteurs. Sans même le programmer, tu peux le faire par Excel, mais je maintiens qu'avec MatLab, c'est beaucoup moins laborieux, pour peu que tu sois familier de la méthode des moindres carrés.
Je vais essayer de te faire un exemple sur Excel avec les données que tu as envoyées.
Par contre, ton histoire de coefficient de friction qui varie, ça m'intéresse. Je fais des recherches en topo de précision et j'étudie de nouveaux systèmes de mesures dont un projet de système de mesure de différences d'altitude par un laché de bille. Je mesure le temps nécessaire au déplacement permettant, par une intégration et connaissant les forces en jeu, de déduire la dénivellée. C'est dans le cas d'un transfert de géométrie entre la surface et le fond d'un tunnel. Aurais-tu des références sur ce point ?
us_30
Messages postés2065Date d'inscriptionlundi 11 avril 2005StatutMembreDernière intervention14 mars 201610 7 juil. 2009 à 23:15
Bonsoir,
Bien compliqué tout cela...
Perso, j'arrive à :
y = 0,00023124*x^2 + 2,12946809*x + 342,20515405
pour la partie courbe, et à :
y=16286,1363636364
pour l'horizontale...
Reste à calculer le point de rencontre... pas très compliqué...
Je ne sais pas si c'est la réponse attendue, ou si tu veux qu'on reprenne un programme qui ne convient pas... mais dans ce cas, je dirais qu'il faut demander à l'auteur...
Ah, oui... C'est une utilisation de base d'Excel...
Nukix
Messages postés5Date d'inscriptionjeudi 2 juillet 2009StatutMembreDernière intervention15 juillet 2009 10 juil. 2009 à 11:02
Salut Thomas,
tu as bien compris pour ce qui est des ballons, en effet j'étudie leur déplacement.
En fait, j'ai déja les position en fonction du temps, je veux juste lisser un peu ça sans perdre les déplacemnts incertains du ballons, chose que j'étudie à vrai dire...
Pourquoi?
car je cherche à optimiser ces courbes pour calculer les meilleurs Cx(coeff de friccion de l'air sur le ballon) et Reynodls, nombre adimensionnel qui dépend du milieu et de la vitesse de déplacement du ballon dans celui ci. Ces renseignements me servent à étudier la mise au plafond des ballons qui parfois dépassent leur plafond car ils arrivent trop vite et éclatent... (le Cx diminue trop avant la mise au plafond)
La majeur partie du temps lorsqu'ils dépassent leur plafond, ils redescendent car l'enveloppe se comprime trop et fait redescendre le ballon qui va osciller autour de l'altitude du plafond pour s'y stabiliser.
C'est pas si facile à expliquer ^^
Pour le moment j'ai des résultats corrects. Je cherche à sortir le polynome calculé.
Le créateur du programme m'a conseillé de je cite :
Pour ta deuxième question, dans la sub main, tu peux
récupérer les coefficients du polynôme, après la
ligne
Call
regression_polynomiale(6, x(), y(), ai_p())
Les coefficients sont les valeurs du vecteur
ai_p()
Pour les récupérer, il suffit par exemple de remplacer
la dernière boucle par celle-ci-dessous :
For j = 1 To
(UBound(ai_p()))
Print #1, Format(ai_p(j), "####0.000")
Next
i
Je vais essayer de faire ce qu'il dit mais je ne comprend pas trop ce que Print #1, Format(ai_p(j), "####0.000") veut dire ^^
Merci de l'interêt que vous portez à mon sujet
Florent
Vous n’avez pas trouvé la réponse que vous recherchez ?
marinmarais
Messages postés104Date d'inscriptionlundi 11 avril 2005StatutMembreDernière intervention16 juillet 20101 13 juil. 2009 à 15:58
Bonjour Florent,
J'ai jeté un coup d'oeil à tes données. Vu la rupture nette à 16100m, une régression polynomiale va te donner des erreurs importantes à cet endroit. Si l'altitude en fonction du temps avait été à peu près parabolique, ça aurait pu jouer mais là...
Bref, j'ai fait une FFT vite fait sur MatLab. Ca marche pas mal... Mon erreur max après modélisation est de 266m (écart-type sur les résidus de 82m). Bah ouais, la FFT ça marche bien sauf aux extrémités.
Si ça te tente, je peux t'envoyer tout ça par mail (fichier excel de résultat et le code MatLab de mon traitement. Dans ce cas, donne moi ton adresse mail par mp.
hogsvoor
Messages postés2Date d'inscriptionlundi 10 novembre 2003StatutMembreDernière intervention16 juillet 2009 13 juil. 2009 à 23:02
bonjour ,
tous roule impec pour moi avec ce programme.
le seul point un peu zarbi est l'appel a des valeurs dans un fichier texte dans la procédure de lancement main
ci dessous comment j'ai modifé la procédure main pour
-charger les données a étudier dans une feuille excel
-choisir les degrés du polynome de regression
Sub main()
Dim degpol As Integer
degpol = 2 ' degrés demandé du polynome de régression , ici 2 ( si un on aura la une droite )
'Appel de la procédure de régression polynômiale :
'le premier argument est le degré du polynôme de tendance choisi, x() et y() sont les données à traiter
'ai_p() est un argument de retour renvoyant les coef ai du polynôme
nbp = 25 ' nombre de points couple x y
ReDim x(nbp)
ReDim y(nbp)
ReDim ai_p(degpol)
For i = 1 To nbp
x(i) = Sheets("prob").Cells(i, 1) ' charge les données x dans une feuille appelée prob de la ligne 1 a 25 colonne 1 à adapter
y(i) = Sheets("prob").Cells(i, 2) ' charge les données y dans une feuille appelée prob de la ligne 1 a 25 colonne 2 à adapter
Next
Call regression_polynomiale(degpol, x(), y(), ai_p()) ' passe les x , y et l'ordre du polynome au programme de calcul
' au retour les coefficients de regression sont le tableau ai_p() de taille degpol
End Sub
sans trop de difficulté il doit être possibles de transformer ces programmes en fonction
si on pouvait avoir en rab le coeff de coreelation ce serait super
marinmarais
Messages postés104Date d'inscriptionlundi 11 avril 2005StatutMembreDernière intervention16 juillet 20101 15 juil. 2009 à 10:58
Salut !
En utilisant la FFT dans ton cas, il ne s'agit pas d'obtenir un spectre en fréquence. Fourier a inventé la transformée juste pour modéliser des fonctions, afin qu'elles soient facilement dérivables pour résoudre des équa diff...
Ici, le but est le même. Tu as l'altitude de ton ballon en fonction du temps. Tu cherches une fonction (pas trop compliquée si possible) modélisant au mieux le mouvement dans le temps de ton ballon.
Et bien tu obtiens comme modèle une somme de cosinus...
Le calcul que j'ai fait, sur MatLab, se déroule en 2 étapes.
- Je ré-échantillonne pour avoir un pas régulier de 60 s et un nombre pair de données
- Je lance la FFT sur les données ré-échantillonnées
Voila le code MatLab. Je suis pas un pro de ce langage, ça peut sans doute être amélioré, mais ça marche !
il n'y pas de limite au nombre de couples de données ( 32000 dans excel 2003 et infini dans excel 2007 sauf la longueur des calculs
ni au nombre de degrés du polynome , mais vu ma (très) longue expérience de modélisateur du vivant , quand un phénomène de régression a deux variables n'est pas toujours pas correctement modélisable par un polynome qui dépasse le 4 ième degré , il y a bien peu de chance que les deux variables soient corrélées ( les mouvements d'un couple de danseur de tangou tango font exception à cette règle)
Dans vba le ficher macro lancé par calcul est la procédure calcul()
j'ai essayé de documenter le truc
précision : n'ayant jamais eu le moindre cours d'excel et de visual basic j'utilise mon propre patois de programmation, désolé pour les pro et les puriste si j'ai un accent et une syntaxe pourris.
mais comme disait le vieux Deng Xiao Ping en 1978 , peu importe la couleur du chat l'important c'est qu'il chope les souris.
nb nikix t'a trouvé où les programmes de calcul parce que celui ou celle qui fait çà , il a fait l'essentiel du travail ce serait bien de lui faire savoir que ça sert à d'autres.
alusipdranreb
Messages postés2Date d'inscriptionjeudi 22 avril 2010StatutMembreDernière intervention25 mai 2010 25 mai 2010 à 09:09
Bonjour.
Je fais des maths sur les nuages de points, dont les nuages plans comme le tien.
Il faut que la fonction numérique f(x) qui représente ton nuage ait un sens physique. Il faut donc lui imposer des contraintes, par exemple sur les coordonnées d'au moins un point du nuage.
Ex:f(x)=7.4447664402*10^-15*x^5-0.000000000149*x^4+0.000000926419*x^3-0.00203323893*x^2+4.1391583705*x-39.4144584247 ne convient pas car pour x=0, ton ballon est enterré à une profondeur de 39.4144584247 m.(x est le temps).Somme des résidus yi-f(xi)=14986.23..
On peut choisir f(x)avec ln.exp.Arctan. ect...
Je te propose un autre polynôme qui satisfait aux contraintes des moindres carrés avec contrainte sur trois points du nuage:(0.200) (4860.16180) (6480.16299)ici.
f(x)=3.42288845778*10^-21*x^7-6.863274651*10^-17*x^6+5.25877769937*10^-13-0.000000001966*x^4+0.000003784328*x^3-0.003400684559*x^2+3.63474227891*x+200. Somme des résidus yi-f(xi)=15834.42.
J'ai pris trois points caractéristiques de ton nuage, dans cet exemple. Si tu as d'autres nuages en dim 2, dim 3..tu peux me les montrer si je peux te citer pour publication dans mon travail. Amicalement. Bernard.