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