Problème d'optimisation avec des données expérimentales

Bonjour,

Je suis étudiante et j'ai effectué des mesures expérimentales et j'aimerais les comparer avec un modèle.
J'ai donc mes données sous forme de liste et le modèle possède 3 inconnues (dont j'ai une mince idée), je recherche la valeur de ces inconnues qui "collerait" le plus avec mes données expérimentales.

Après quelques recherches internet, j'ai découvert le module scipy.optimize mais j'ai du mal à comprendre les arguments qu'il faut y mettre et à faire mon choix parmi minimize / fmin et compagnie... Je ne comprends pas bien comment obtenir ce que je souhaite avec ces fonctions.

Si vous avez la moindre idée, n'hésitez pas!

Bonjour

Il faudrait déjà que tu définisses ce que tu entends par "collerait". Perso j'arrive pas à visualiser comment faire coller "une" inconnue avec une "liste" de valeurs (ou plutôt moi je visualise bien quelque chose mais rien ne garantit que ma vision de l'opération corresponde à la tienne). C'est d'ailleurs sûrement la même chose avec "optimize". Peut-être que cette fonction ne correspond pas non plus à ta vision de "collerait"...

Salut,

Le choix de la méthode d'optimisation va dépendre du type de loi que tu cherche à "coller" (linéaire, non linéaire?) et des contraintes eventuelles que tu as (les variables sont-elles bornées?).

Est-ce que tu peux nous donner cette loi? On pourra t'orienter vers une méthode en particulier et t'expliquer les paramètres au besoin.

Le principe générale reste tout de même identique: on define un fonction "objectif" qui doit être le plus proche possible de zero (ici on a quelque chose comme la somme des écarts au carré entre les données expérimentales et les données provenant de la loi le plus proche possible de zero).

ça devrait le faire.

J

Pour être plus précise, j'ai des mesures de masses faites à différents instants, j'ai donc une liste avec les valeurs de masse en teneur en eau et une liste avec les valeurs de temps. J'ai tracé l'évolution de la masse en fonction du temps et je cherche à faire une modélisation de la loi de Fick 1D.
La loi de Fick 1D permet d'obtenir la masse en eau en fonction du temps : M(t)=Msat*4/e*sqrt(D*t/pi) avec e l'épaisseur de l'échantillon (connue)
Le but serait donc de trouver les valeurs de Msat (Masse à saturation) et D (coefficient de diffusion) pour que les valeurs de M(t) soient le plus proche possible des valeurs expérimentales de masse.

Je vous donne un bout de programme calculant une différence entre la fonction théorique et les valeurs expérimentales. Cependant, il ne marche pas, je pense que je n'ai pas bien compris le fonctionnement de optimize.

teneur_eau1 et temps sont des listes de même longueur

Code:

---

1 2 3 4 5 6 7

def f(D,Msat): res=0 for i in range(0,len(temps)): res+=pow(teneur_eau1[i]-Msat/e*pow(D*temps[i]/m.pi,0.5),2) return(res) xo=np.array([7.89495962707063*1E-005,teneur_eau1[-1]]) resultat=opt.minimize(f,xo,(D,Msat),method='nelder-mead')

---

Hmm. Cela risqué d'être compliqué. La formulation est non seulement non-linéaire mais en plus elle n'est pas assez "contrainte". Trouver un coefficient A = Msat*racine(D) c'est assez facile, mais Msat et D peuvent prendre n'importe quelles valeurs...

Voici un exemple utilisant la méthode "SLSQP" pour Sequential Least SQuares Programming.

On crée une function fick1D:

Code:

---

1 2 3

def fick1D(saturated_mass, diff_coeff, time, thickness): """Diffusion funtion (1D)""" return saturated_mass*4/thickness*(diff_coeff*time/np.pi)**0.5

---

On definit une autre function en charge de qualifier notre fitting. Ici on calcul la somme des écart au carré. C'est le résultat de cette fonction que l'on cherche à minimiser:

Code:

---

1 2 3 4 5

def residuals(x, measured_mass, measured_time, thickness): """Residual function. Return sum of square residuals""" saturated_mass, diff_coeff = x fitted_mass = fick1D(saturated_mass, diff_coeff, measured_time, thickness) return np.sum((measured_mass - fitted_mass)**2)

---

Afin d'illustrer un cas concret je génère un jeu de données plus ou moins réaliste...

Code:

---

1 2 3 4 5 6 7 8

# SOME DATA FOR TESTING PURPOSE TRUE_DIFF_COEFF = 1.0E-04 SATURATED_MASS = 20.0 THICKNESS = 1.0   t = np.arange(0, 10800, 1) true_mass = fick1D(SATURATED_MASS, TRUE_DIFF_COEFF, t, THICKNESS) measured_mass = np.random.normal(0, 0.5, len(t)) + true_mass

---

Note que j'ai créé un série fictive de valeurs mesurées en ajoutant de l'aléatoire. Maintenant l'optimisation:

Code:

---

1 2 3 4 5 6

# Optimization x0 = (1.0, 2.5E-04) # (saturated_mass, diff_coeff) args = (measured_mass, t, THICKNESS) res = optimize.minimize(residuals, x0, args=args, method='SLSQP') fitted_saturated_mass, fitted_diff_coeff = res.x fitted_mass = fick1D(fitted_saturated_mass, fitted_diff_coeff, t, THICKNESS)

---

Ici x c'est ce sur quoi on peut jouer pour optimiser. C'est (saturated_mass, diff_coeff). J'ajoute des arguments à la function d'optimisation, ces arguments sont utilisés par la function residuals, ce sont ici la time series de la mesure.

On peut verifier que le résultat est bon et mauvais:

Code:

---

1 2 3 4 5 6 7 8 9 10 11 12

fig = plt.figure() ax = fig.add_subplot(111) ax.plot(t, true_mass, label='True Mass', color='black') ax.scatter(t, measured_mass, label='Measured Mass', color='red') ax.plot(t, fitted_mass, label='Fitted Mass', color='blue', linewidth=2) ax.set_xlabel('TIME (S)') ax.set_ylabel('MASS (KG)') ax.set_title('SATURATED MASS = %0.3e - DIFFUSION COEFF = %0.3e' % \ (fitted_saturated_mass, fitted_diff_coeff)) ax.legend(loc='upper right', prop={'size':11,}) plt.grid() plt.show()

---

Le graph ressemble à ceci:
Pièce jointe 178172
Tu peux voir que bien que l'on fitte très bien.Les valeurs de msat et D sont dans les choux....


Ju

Merci beaucoup pour cette réponse très complète!

J'avais cru comprendre qu'on pouvait borner les valeurs des variables recherchées (ici Msat et D) grâce aux options de la fonction optimize.minimize, notamment "bounds". Cela ne permettrait pas de résoudre le problème dans ce cas là?
J'ai lu dans de nombreuses thèses sur le phénomène de diffusion que la méthode la plus fréquemment utilisée pour résoudre ce type d'équation était nelder-mead. Y-a-il une grande différence? L'équation que j'ai donné pour la résolution de Fick 1D est simplifié, on a normalement des sommes infinies, peut-être que cela fausse un peu les calculs...

J'ai fait tourner l'algorithme en modifiant quelques fonctions et j'obtiens un résultat "NaN" en mode SLSQP et un résultat réel mais incohérent avec la méthode nelder-mead ..

Pièce jointe 178199

La courbe théorique modélise très mal les points expérimentaux...

Salut,

Citation:

---

Envoyé par demaina

J'avais cru comprendre qu'on pouvait borner les valeurs des variables recherchées (ici Msat et D) grâce aux options de la fonction optimize.minimize, notamment "bounds". Cela ne permettrait pas de résoudre le problème dans ce cas là?

---

Effectivement ça peut aider. J'avoue avoir essayé vite fait hier avec

Code:

---

1 2 3

bnds = ((0.0, 1000.0), (1.0E-04, 5.0E-04)) res = optimize.minimize(residuals, x0, args=args, method='SLSQP', bounds=bnds)

---

Mais j'ai le message d'erreur suivant:

Code:

---

'Inequality constraints incompatible'

---

Je n'ai peut-être pas les yeux en face des trous mais je ne vois pas en quoi mes contraintes ne sont pas bonnes...

Sinon, tes mesures semblent bien suivre une droite en function de racine(t), mais la courbe fittée elle semble être trace en function de t seulement. Je n'ai pas d'avis sur la method d'optimization.. Il faudrait soit reposter un message dans une section math ou demander à déplacer se sujet dans la section "scientifique".

J

En effet, ma liste de temps était en racine seconde et non en seconde ce qui faussait tous les résultats... Et ça fait maintenant quelques semaines que je cherchais d'où provenait cette erreur! Voici le résultat obtenue avec cette simple modification :

Pièce jointe 178256

Si quelqu'un a une solution pour mettre les graduations des abscisses en notation scientifique, ça m'intéresserait!

Certes, la Masse à saturation pas forcément très cohérente, mais la courbe théorique modélise très bien les points expérimentaux! Je vais essayer de rajouter des bornes pour les variables pour voir si je peux obtenir une masse un peu plus cohérente ou tester un modèle moins simplificateur avec des sommes infinies. Un grand merci à toi pour toutes ces explications sur les optimisations et pour m'avoir grandement aidé!

Pour passer en notation scientifique, tu peux essayer ceci:

Code:

---

plt.ticklabel_format(axis='x', style='sci', scilimits=(1,4))

---

Comme je le disais plus haut, le problème c'est qu'il manque une contrainte pour obtenir un bon ordre de grandeur des variables Msat et D. Actuellement il y a une infinité de solutions... Même avec des contraintes (bounds) cela ne suffira pas. Il faut peut-être supposer une donnée comme connue, par exemple le coefficient de diffusion (si c'est de l'eau, la valeur est connue). Après je ne connais pas la problèmatique sous-jacente et tu peux ne pas pouvoir figer une variable. Ou alors tu peux utiliser des données issue d'un autre système de mesure. Comme un thermomètre, ou autre, pour estimer D...

J