Résolution d'équations non linéaires en python

Il y a deux façons de le faire.

1. utiliser un solveur non linéaire
2. Linéarisez le problème et résolvez - le dans le sens des moindres carrés

Setup

donc, comme je comprends votre question, vous savez F, a, b, et c à 4 points différents, et vous voulez inverser pour les paramètres du modèle X, Y, et Z. Nous avons 3 inconnues et 4 points de données observées, donc le problème est surdéterminé. Par conséquent, nous résoudrons dans les moindres carrés sens.

il est plus courant d'utiliser la terminologie opposée dans ce cas, alors retournons votre équation. Au lieu de:

F_i = X^2 + a_i Y^2 + b_i X Y cosZ + c_i X Y sinZ

nous allons écrire:

F_i = a^2 + X_i b^2 + Y_i a b cos(c) + Z_i a b sin(c)

Où nous savons F,X,Y et Z en 4 points différents (par exemple F_0, F_1, ... F_i).

nous changeons simplement le nom des variables, pas l'équation elle-même. (Ceci est plus pour ma facilité de penser que n'importe quoi d'autre.)

Linéaire Solution

Il est possible de linéariser cette équation. Vous pouvez facilement résoudre pour a^2,b^2,a b cos(c) et a b sin(c). Pour rendre cela un peu plus facile, ré-étiquetons les choses une fois de plus:

d = a^2
e = b^2
f = a b cos(c)
g = a b sin(c)

Maintenant, l'équation est beaucoup plus simple: F_i = d + e X_i + f Y_i + g Z_i. Il est facile de faire une inversion linéaire des moindres carrés pour d,e,f et g. Nous pouvons alors obtenir a,b et c:

a = sqrt(d)
b = sqrt(e)
c = arctan(g/f)

OK, écrivons ceci en forme de matrice. Nous allons traduire 4 observations de (le code que nous allons écrire va prendre n'importe quel nombre d'observations, mais gardons le concret pour le moment):

F_i = d + e X_i + f Y_i + g Z_i

En:

|F_0|   |1, X_0, Y_0, Z_0|   |d|
|F_1| = |1, X_1, Y_1, Z_1| * |e|
|F_2|   |1, X_2, Y_2, Z_2|   |f|
|F_3|   |1, X_3, Y_3, Z_3|   |g|

Ou: F = G * m (je suis un geophysist, de sorte que nous utilisons G pour "fonctions de Green" et m pour "Paramètres du modèle". Généralement, on utilisera d pour "data" au lieu de F.)

en python, cela se traduirait à:

def invert(f, x, y, z):
    G = np.vstack([np.ones_like(x), x, y, z]).T
    m, _, _, _ = np.linalg.lstsq(G, f)

    d, e, f, g = m
    a = np.sqrt(d)
    b = np.sqrt(e)
    c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
    return a, b, c

Non-linéaire de la Solution

vous pouvez aussi résoudre ce problème en utilisant scipy.optimize, comme @Joe l'a suggéré. La fonction la plus accessible en scipy.optimizescipy.optimize.curve_fit qui utilise une méthode Levenberg-Marquardt par défaut.

Levenberg-Marquardt est un algorithme de" hill climbing " (bien, il descend en pente, dans ce cas, mais le terme est utilisé de toute façon). Dans un sens, vous faites une première estimation des paramètres du modèle (tous, par défaut dans scipy.optimize) et suivre la la pente de observed - predicted dans votre espace de paramètre descente vers le bas.

Avertissement: choisir la bonne méthode d'inversion non linéaire, deviner au départ, et ajuster les paramètres de la méthode est un "art sombre". Vous ne l'apprenez qu'en le faisant, et il y a beaucoup de situations où les choses ne fonctionneront pas correctement. Levenberg-Marquardt est une bonne méthode générale si votre espace de paramètre est assez lisse (celui-ci devrait l'être). Il y a beaucoup d'autres (y compris la génétique des algorithmes, réseaux neuronaux, etc. en plus de méthodes plus courantes comme le recuit simulé) qui sont meilleurs dans d'autres situations. Je ne vais pas plonger dans la partie ici.

il y a un gotcha commun que certains outils d'optimisation tentent de corriger pour cela scipy.optimize n'essaie pas de gérer. Si les paramètres de votre modèle ont des valeurs différentes (par exemple a=1, b=1000, c=1e-8), vous aurez besoin de réviser les choses pour qu'elles soient de la même ampleur. Sinon,scipy.optimize'"escalade" des algorithmes (comme LM) ne calculera pas avec précision l'estimation du gradient local, et donnera des résultats extrêmement inexacts. Pour l'instant, je suis en supposant que a,b et c ont des grandeurs relativement semblables. En outre, soyez conscient que essentiellement toutes les méthodes non-linéaires vous exigent de faire une première supposition, et sont sensibles à cette supposition. Je le laisse en dessous (passez-le comme le p0 kwarg curve_fit) parce que la valeur par défaut a, b, c = 1, 1, 1 est une estimation assez exacte pour a, b, c = 3, 2, 1.

Avec les mises en garde de la route, curve_fit s'attend à passer une fonction, un ensemble de points où les observations ont été faites (en un seul ndim x npoints tableau) et les valeurs observées.

Donc, si nous écrivons la fonction comme ceci:

def func(x, y, z, a, b, c):
    f = (a**2
         + x * b**2
         + y * a * b * np.cos(c)
         + z * a * b * np.sin(c))
    return f

nous aurons besoin de l'envelopper pour accepter des arguments légèrement différents avant de le passer à curve_fit.

En bref:

def nonlinear_invert(f, x, y, z):
    def wrapped_func(observation_points, a, b, c):
        x, y, z = observation_points
        return func(x, y, z, a, b, c)

    xdata = np.vstack([x, y, z])
    model, cov = opt.curve_fit(wrapped_func, xdata, f)
    return model

exemple autonome des deux méthodes:

pour vous donner une implémentation complète, voici un exemple qui

1. génère des points distribués au hasard pour évaluer la fonction sur,
2. évalue la fonction sur ces points (en utilisant les paramètres du modèle set),
3. ajoute du bruit aux résultats,
4. et ensuite s'inverse pour les paramètres du modèle en utilisant à la fois les méthodes linéaires et non linéaires décrites ci-dessus.

---

import numpy as np
import scipy.optimize as opt

def main():
    nobservations = 4
    a, b, c = 3.0, 2.0, 1.0
    f, x, y, z = generate_data(nobservations, a, b, c)

    print 'Linear results (should be {}, {}, {}):'.format(a, b, c)
    print linear_invert(f, x, y, z)

    print 'Non-linear results (should be {}, {}, {}):'.format(a, b, c)
    print nonlinear_invert(f, x, y, z)

def generate_data(nobservations, a, b, c, noise_level=0.01):
    x, y, z = np.random.random((3, nobservations))
    noise = noise_level * np.random.normal(0, noise_level, nobservations)
    f = func(x, y, z, a, b, c) + noise
    return f, x, y, z

def func(x, y, z, a, b, c):
    f = (a**2
         + x * b**2
         + y * a * b * np.cos(c)
         + z * a * b * np.sin(c))
    return f

def linear_invert(f, x, y, z):
    G = np.vstack([np.ones_like(x), x, y, z]).T
    m, _, _, _ = np.linalg.lstsq(G, f)

    d, e, f, g = m
    a = np.sqrt(d)
    b = np.sqrt(e)
    c = np.arctan2(g, f) # Note that `c` will be in radians, not degrees
    return a, b, c

def nonlinear_invert(f, x, y, z):
    # "curve_fit" expects the function to take a slightly different form...
    def wrapped_func(observation_points, a, b, c):
        x, y, z = observation_points
        return func(x, y, z, a, b, c)

    xdata = np.vstack([x, y, z])
    model, cov = opt.curve_fit(wrapped_func, xdata, f)
    return model

main()