Discussion :

[demaina]

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!

[Sve@r]

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"...

[Julien N]

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

[demaina]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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')

[Julien N]

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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:

Tu peux voir que bien que l'on fitte très bien.Les valeurs de msat et D sont dans les choux....


Ju

[demaina]

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 ..



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

[Julien N]

Salut,

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

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 :

Sélectionner tout - Visualiser dans une fenêtre à part

'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

[demaina]

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 :



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é!

[Julien N]

Pour passer en notation scientifique, tu peux essayer ceci:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

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